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EXECUTIVE  SUMMARY  -  PHASE  I 


In  the  design  and  operation  of  large  space  structures,  the  effects 
of  time  dependent  and  random  vibrations,  heating  and  cooling,  control 
forces,  etc.,  lead  to  severe  analytical  problems  because  of  nonlinearities, 
time-varying  effects  such  as  self-shadowing  In  trusses,  couplings,  etc. 
These  problems  lead  to  model  equations  which  can  be  differential,  partial 
differential,  or  integro-differential  equations  involving  complex 
nonlinearities  and  stochastic  parameters. 

Our  objective  is  to  identify  these  problems  and  to  demonstrate  that 
the  Adomian  decomposition  method  provides  a  new  useful  approach  to  such 
problems  and  that  the  approach  provides  physically  more  realistic  and 
correct  solutions,  and,  further,  that  tremendous  computational  effort 
can  be  avoided. 

The  saving  of  computer  time  results  from  the  avoidance  of  discretized 
methods  and  the  ability  to  obtain  a  continuous  analytical  solution.  This 
not  only  saves  computation  time  but  makes  it  easier  to  see  functional 
dependences. 

The  fact  that  solutions  are  necessarily  more  correct  arises  from 
avoidance  of  the  conventional  linearization  and  restrictive  assumptions 
normally  required.  The  ease  and  accuracy  of  the  method  is  demonstrated  in 
a  variety  of  selected  examples  and  problems. 

In  addition  to  the  fulfillment  of  the  above  contract  objective, 
some  very  significant  other  and  additional  new  insights  and  progress  have 
been  achieved  to  be  expanded  on  in  Phase  II.  We  have  established  that 
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I  parabolic,  elliptic,  and  hyperbolic  equations  can  all  be  solved  by  the 

I  decomposition  method  (there  Is  no  other  method  that  can  do  this 

K  analytically. 

y  Secondly  we  have  seen  as  a  direct  result  of  this  study  that  we  can 

r  generalize  Kalman  filters  to  a  continuous  nonlinear  stochastic  case 

I  without  white  noise  limitations  and  that  we  can  generalize  control  theory 

[  analytically  to  nonlinear  stochastic  multidimensional  control  theory  for 

p  space  structures,  an  insight  of  the  greatest  significance  to  be  an 

Important  area  of  our  Phase  II  proposed  work. 
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1)  In  the  Journal  of  American  Institute  of  Aeronautics  and  Astronautics 
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in  1988. 

3)  In  the  Journal  oj  Mathematical  Analyses  and  Applications 

"A  Review  of  the  Decomposition  Method,"  G.  Adomian,  to 
appear  Sept.  1988. 
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Principal  accomplishments: 

1)  The  contract  statement  of  work  has  been  fulfilled  in  establishing 
that  the  decomposition  method  is  valuable  in  the  basic  types  of 
typical  model  equations  that  will  arise  in  analysis  of  large 
space  structures  such  as  vibration,  heating,  etc.  It  has  been 
shown  that  the  decomposition  method  solves  typical  nonlinear  and 
stochastic  and  multidimensional  model  equations  without  customary 
restrictive  and  limiting  assumptions  or  procedures  involving 
excessive  computation. 

2)  Additional  insights  as  a  result  of  this  study  are: 

a)  the  procedure  will  apply  to  parabolic,  elliptic,  and  hyperboli 
equations,  a  unique  feature  and  advantage  of  the  decomposition 
method,  which  means  rather  global  applicability  to  aerospace 
problems. 

b)  Control  theory  can  be  generalized  to  nonlinear  stochastic  and 
multidimensional  cases  for  application  to  space  structures 
and  will  be  studied  in  detail  in  Phase  II.  Kalman  filter 
theory  can  similarly  be  generalized  to  the  case  of  nonlinear 
stochastic  operators  in  both  s.ate  and  measurement  equations. 


INTRODUCTION 


^The  large  structures  contemplated  would  be  constructed  in  space. 
Because  of  the  limitations  on  launching  massive  payloads,  it  is  clear 
that  these  structures  will  be  made  of  lightweight  material  and  will 
necessarily  be  flexible  and  easily  excited  into  vibrations. 

Analytical  problems  will  arise  in  designing  large  space  structures 
in  which  physically  realistic  and  accurate  solutions  will  be  critical. 
Such  designs  must  consider  weight,  sizes,  stiffness,  thermal  and 
mechanical  distortions,  stresses  due  to  gravity  and  positioning  thrusts. 
Some  specific  analytical  problems  will  involve  vibration,  heating  and 
cooling,  multidimensional  control,  and  structural  problems  arising  from 
random  support  motion  and  random  fluctuations  of  the  system  dynamic 
properties.  - — & fjr  ?  S)  ^ 

In  large  space  structures,  vibrations  can  result  from  machinery, 
thermal  transients,  and  operational  maneuvers.  The  response  to  such 
vibrations  depends  on  structural  parameters  (stiffness,  mass,  damping) 
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and  the  ex  jrnal  excitations  and  on  the  control  system. 

Vibration  is  hardly  the  only  problem.  There  will  be  effects  of 
severe  heating  (e.g.,  on  sun-illuminated  parts)  and  severe  cooling  on 
dark  sides  with  resulting  exapnsions,  contractions,  and  stresses  on  the 
structure.  Additionally  there  will  be  the  problem  of  nonlinear  stochastic 
multidimensional  control  of  a  large  complex  structure  with  limited 
rigidity,  and  the  interactions  through  the  central  system  as  a  result  of 
frequencies  introduced  in  the  control  loop  through  rotation,  bending, 
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expansion  and  contraction.  The  means  employed  for  passive  damping  also 
complicate  matters  since  modifications  of  mass,  stiffness,  and  damping 
affect  the  control  design. 

If  a  flexible  structure  were  designed  to  be  rotated  at  a  constant 

angular  velocity,  then  complex  Coriolis  and  centripetal  accelerations 

as  well  as  accelerations  due  to  structural  deformations  would  have  to 

be  taken  into  account.  Thus  the  rotational  motions  are  a  factor  in 

dynamic  analysis.  Ordinarily,  such  analyses  immediately  rely  on 

linearizing  the  governing  equations  of  motion.  The  decomposition  method, 

on  the  other  hand,  allows  us  to  avoid  this  restrictive  procedure. 

The  fact  that  the  structural  materials  will  be  nonhomogeneous 

composites  will  result  in  random  fluctuations  in  dynamic  properties. 

Furthermore,  once  the  structure  begins  to  vibrate,  structural  parameters 

may  fluctuate.  Thus  one  must  be  able  to  deal  with  cases  in  which  the 

stiffness  or  damping  are  random  functions  of  time  rather  than  constants. 

The  equations  involved  will  then  be  stochastic  differential  equations  in 

★ 

the  sense  used  by  the  principal  investigator  in  the  referenced  works  , 
i.e.,  differential  equations  with  stochastic  process  coefficients  rather 
than  the  much  simpler  case  with  only  stochastic  inputs  -  also  called  s.d.e. 
by  many  investigators. 


In  the  so-called  "disordered  systems"  where  parameters  are  random 
variables  we  have  random  boundary-value  problems.  The  decomposition 
method  can  deal  with  boundary  operator  equations  which  are  random. 
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nonlinear,  or  even  coupled. 


★ 

See  particularly  references  [1,2]. 

** 

Assuming  that  systems  parameters  are  merely  uncertain  rather  than 
fluctuating  in  time. 
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One  should  consider  various  types  of  nonlinearities  in  order  to  gain 
insights  into  possible  unusual  or  growing  structural  vibrations  that  might 
lead  to  failure.  Thus  our  ability,  when  using  decomposition,  to  consider 
strong  and  composite  nonlinearities  is  of  extreme  importance. 

Introduction  to  the  Decomposition  Method: 

Real  physical  systems  are  nonlinear  and  stochastic;  linearity  and 
determinism  are  special  cases.  Realistic  modeling  will  lead  to  systems 
of  ordinary  or  partial  (nonlinear  stochastic)  differential  equations. 

Our  model  equations  will  involve  nonlinearities  and/or  stochastic 
parameters,  inputs,  and  initial /boundary  conditions. 

The  methods  that  have  been  available  so  far  to  solve  such  systems 
have  limitations  which  require  close  examination.  It  has  become 
standard  to  linearize  nonlinear  equations  or  to  assume  they  are  "weakly 
nonlinear."  In  the  stochastic  case,  resort  to  statistical  linearization 
is  a  common  procedure.  Stochasticity  in  parameters  is  quite  commonly 
ignored.  Inputs  which  we  know  to  be  stochastic  processes  are  assumed 
tc  be  delta-correlated  or  to  represent  "small"  fluctuations.  Thus  unless 
systems  are  essentially  linear  and  deterministic,  i.e.,  unless  nonlinearity 
and  stochasticity  are  not  significant  factors,  what  is  commonly  being 
done  -  whether  apparent  or  not  -  is  to  change  the  system  physical  model 
to  a  mathematically  more  tractable  one  within  the  capability  of  the 
available  mathematics.  The  contrast  between  "closed-form  solutions" 
and  "approximate  solutions"  has  been  exaggerated  and,  sometimes,  not  even 
clearly  understood.  The  reason  is  simple.  All  modelling  is  approximation. 
In  modelling,  we  always  have  a  compromise  between  realistic  modeling  and 


tractability  In  our  available  methodologies.  The  model  actually  solved 


by  conventional  methods  is  often  not  the  model  we  had  wanted  to  solve 


which  was,  of  course,  already  and  necessarily  an  approximation  to  reality. 


With  the  model  forced  into  a  tractable  form  by  restrictive  assumptions 


and  techniques  such  as  linearization,  one  can  then  find  "nice"  closed-form 


solutions.  However  this  only  means  we  have  an  identifiable  or  well-known 


series  such  as  sin  x  =  x  -  xJ/3i  +  •••  as  the  solution  for  a  mathematically 


simplified  problem. 


Our  method  of  attack  is  a  continuous  approximation  method  which  yields 


a  rapidly  converging  series  solution  (which  can  sometimes  be  summed  into 


a  familiar  closed  form). 


What  is  significant  is  that  it  does  not  require  linearization.  It 
solves  nonlinear  problems.  We  do  not  limit  stochastic  cases  to  perturbation 


theory,  Wiener  processes,  delta-correlated  processes,  Markovian  assumptions. 


closure  approximations,  or  quasi -monochromatic  assumptions. 


This  method  is  computationally  convenient.  It  does  not  require 
discretization  into  grids  with  its  resulting  massive  computation  and 
approximation  between  computing  points.  It  yields  a  continuous  analytic 


solution  into  which  we  then  put  numbers.  We  can  then  see  the  rapid 


stabilization  of  the  solution  as  terms  are  calculated. 


Such  a  methodology  -  we  have  called  it  the  "decomposition  method"  - 
which  yields  solution  for  a  model  retaining  actual  nonlinearity  and 
stochasticity,  seems  preferable  to  a  nice  mathematical  solution  -  a  so-called 


exact  solution  -  to  an  artificial  or  simplified  problem  in  which  one 


A  proof  of  the  convergence  has  now  been  made  by  the  principal 
investigator  and  is  in  publication. 
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linearizes  or  assumes  unphysical  processes.  For  example,  most  nonlinear 
stochastic  oscillator  problems  assume  only  stochastic  excitation  and 
require  hierarchy  methods  and  closure  approximations.  However,  such 

** 

assumptions  have  been  shown  to  be  equivalent  to  perturbation  methods. 

In  a  case  where  the  perturbation  method  is  not  adequate  to  the  problem, 
e.g.,  where  fluctuations  are  large,  decomposition  offers  a  more  accurate 
solution.  Decomposition  is  the  more  general  method  always  resulting  in 

•toft 

a  correct  solution.  In  the  cases  where  perturbative  methods  have  been 
adequate,  or  linear  and/or  deterministic  solutions  would  be  sufficient, 
decomposition  will  yield  the  same  solution. 

In  large  (flexible)  space  structures  a  central  problem  will  be 
random,  as  well  as  periodic,  vibrations.  For  random  vibrations,  often 
the  usual  treatments,  such  as  the  Fokker-Planck-Kolmogorov,  perturbation, 
statistical  linearization,  Gram-Charlier  expansions,  averaging  and 
cumulant  closure  methods,  etc.,  are  less  than  desirable  for  many  reasons 
(particularly  more  limitations  and  less  accuracy  in  many  real  problems). 
Restrictive  approximations  that  are  made,  e.g.,  white  noise  excitation 
and  deterministic  parameters,  Wiener  processes,  Gaussian  behavior,  etc., 
may  well  result  in  mathematical  solutions  which  are  incomplete  or  which 
deviate  from  the  actual  physical  behavior. 

We  believe  the  decomposition  method  offers  a  new  approach  with  a 
significant  potential  for  major  contributions  not  possible  otherwise 
because  of  the  avoidance  of  restrictions  necessary  in  other  approaches. 


* 

See  reference  [1]. 

** 

See  references  [1]  and  [2]  . 
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ILLUSTRATIONS  OF  THE  DECOMPOSITION  METHOD 


Solution  of  nonlinear  equations  arising  in  the  modeling  of  a  physical 
system  generally  begins  with  some  form  of  linearization,  assumptions  of 
weak  nonlinearity,  and  "smallness."  Yet  physical  systems  are  nonlinear, 
and  real  systems  generally  involve  random  fluctuations.  The  general  case 
is  a  nonlinear  stochastic  system,  and  "linear"  and  "deterministic"  are 
special  cases.  We  aim  to  solve  physical  problems  realistically,  not 
simply  artificial  versions  of  those  problems  selected  so  as  to  make  the 
mathematics  tractable.  It  is  true  that  we  have  gone  a  long  way  with 
the  earlier  methods,  and  in  many  problems  they  are  completely  adequate; 
in  others,  they  are  a  good  first  approximation  to  gain  some  insight.  Yet, 
one  can  find  cases  in  which  they  are  not  adequate;  clearly  not  all  systems 
are  linear  and  deterministic.  Further,  all  effects  and  behaviors  in 
complex  systems  are  not  instantaneous,  and  a  serious  attempt  to  explain 
behavior  of  a  complex  dynamical  system,  and  eventually  to  control  it, 
must  sometimes  consider  delays  or  retarded  effects. 

To  solve  such  problems,  we  must  be  able  to  solve  equations  or  systems 
of  equations  which  may  involve  differential  or  partial  differential 
operators,  linear  or  nonlinear,  deterministic  or  stochastic,  and  possibly 
involving  delays. 

As  ambitious  as  it  may  appear,  it  is  fortuitous  that  the  decomposition 
method  appears  to  be  capable  of  such  solutions  in  a  fairly  wide  class  of 
problems.  The  superficial  resemblance  of  this  method  to  s:me  other 
methods  can  be  misleading;  the  decomposition  method  solves  problems  not 
solvable  by  other  methods,  or  which  are  only  solvable  with  much  more 


difficulty  or  computation.  The  method  is  an  "approximation"  method,  not 
a  "closed  form"  solution.  The  usual  significance  of  these  terms  is  that 
in  the  one  case,  we  have  an  exact  answer  and  in  the  other  an  approximate 
one.  Clearly,  however,  a  method  of  solution  which  changes  the  problem  to 
a  different ,  easier  mathematical  problem  and  then  solves  it  exactly  is 
not  to  be  preferred  to  one  In  which  the  actual  nonlinear  and/or  stochastic 
model  is  treated  with  an  "approximate"  method  which  provides  accurate, 
rapidly  convergent,  and  computable  series  of  terms. 

Let  us  consider  how  the  decomposition  method  may  change  the  situation. 
Consider  first  a  deterministic  ordinary  differential  equation  Fu  =  x(t) 
where  Fu  is  given  as  the  sum  of  a  linear  term  and  a  nonlinear  term. 

The  linear  term  is  Lu  +  Ru,  where  L  is  the  highest-ordered  derivative 
and  R  is  the  remainder  of  the  linear  operator.  If  Nu  is  the  nonlinear 
term,  Lu  +  Ru  +  Nu  »  x(t),  which  we  write  Lu  =  x  -  Ru  -  Nu.  If  this 
is  an  initial-condition  problem,  then  for  L  =  dn/dtn,  we  define  L-1 
as  the  n-fold  definite  integration  from  0  to  t. 

From  L_1Lu  =  L-1x  -  L-1Ru  -  L”^Nu  we  obtain 

u  =  uQ  -  L_1Ru  -  L_1Nu, 

where  Uq  =  L~\  +  u(0)  +  tu'(O)  +  •••  +  tn"^u^n~^  (0)/(n-l ) !  Now  u  is 

OO  00 

decomposed  into  \  un  and  Nu  is  written  \  An(ug,u^ . un) 

where  the  An  are  the  Adomian  polynomials  generated  fro  the  specific 
nonlinearity  to  represent  it  exactly  [23].  These  polynomials  have  now 
been  discussed  in  numerous  papers  and  books  and  are  listed  for  convenience 


at  the  end  of  this  section  .  Now,  If  Nu  =  f(u)  is  an  analytic  function, 

oo 

we  can  write  f(u)  =  T  An  in  a  series  which  generally  converges  quite 

n=0  n 


ine 


**  -l  ” 

rapidly  .  Now  u  =  un  -  LR  V  u„  -  L"  T  A„  and  we  determii 

0  n=0  n  n=0  n 

un+l  =  'L  1  Rup  -  L  ^An  for  n  ^  0  to  find  all  components  of  u,  provided 

we  have  the  appropriate  integrabil ity.  The  An  depend  on  the  specific 

nonlinearity.  Further,  Ag  depends  only  on  Ug;  A^  depends  on  Ug,  u^ , 

etc.  This  makes  the  solution  calculable. 

If  boundary  conditions  are  specified  instead  of  initial  conditions, 

we  would  allow  L  ^  to  be  n-fold  indefinite  integrations,  then  evaluate 

the  constants. 

The  existence,  uniqueness,  and  other  properties  of  the  solution  depend 
on  the  terms  in  uQ.  In  the  linear  case  where  Nu  vanishes  we  have 

u  =  Ug  -  L_1Rug  -  L'1Ru1  -  •••  =  Ug  -  L_1RUg  +  (L_1R) (L_1R)Ug  -  •••  or 

oo 

u  =  I  ( -1 )n( L  'R)nun.  If  the  given  conditions  are  zero,  un  =  L_1g  and 
n=0  u  u 

oo 

u  =  £  (-1 )n(L  ^ R) nL  ^g.  Thus  Fu  =  g  becomes  u  =  F"^g  where  the 

n=0 

-CO 

inverse  is  F  =  ][  (-1  )n(L~ *R)nL-^ .  This  is,  of  course,  a  succinct 
n=0 

introduction;  general  cases  are  discussed  elsewhere  [1-3]. 


* 

We  have  listed  A  for  n  =  0,1 ,2, . . . ,10.  References  [1,2,3,23] 
show  derivations.  n 
** 

The  1  An  Is  a  rearranged  generalized  Taylor  series.  See  Nonlinear 
Stochastic  Systems  Theory  and  Applications  to  Physics,  G.  Adomian,  Reidel 


publishing  Co.,  1988. 


Example:  As  a  simple  (linear  deterministic)  example,  consider  the 

well-known  Airy's  Equation  y"  -  ty  *  0  which  arise  in  applications  in 

mathematical  physics  and  are  useful  for  asymptotic  representations  of  various 

special  functions.  Assume  conditions  y (0)  =  1  and  y'(0)  =  1  and  write 

2  2 

the  equation  in  the  form  Ly  -  Ry  =  0  with  L  -  d  /dt  ,  R  =  t. 


L"1  =  [  (f  [*]dt)dt.  Then  operating  with  L‘ 

Jn 


we  obtain 


'0  >s> 


y(t)  =  y(0)  +  ty'(0)  +  L~  Ry.  Thus 


yn  =  1  +  t 


,-lDu  t3  .  t4  _  l*t3  2*t4 

L  Ry0  -  L  t(l  +  t)  -  273  +  -  -JT-  +  -jr- 


t6  t7  _  l*4*6*t6  J  2*5*t7 

y2  =  L  ^1  =  - TT 


1*4*7. ..(3n-2)t3n  .  2*5*8... (3n-l)t3n+1 
- (3n71  - + - PStT 


The  sum  7  y  is  the  solution  y.  For  t  *  0.1,  one  term  is  sufficient 
n=0  n 

for  six  decimal  accuracy.  For  t  =  1.0,  three  terms  are  sufficient  to  at 
least  three  digits. 


Example:  Consider  the  equation  d^u/dx^  -  kxpu  =  g  with  u(l)  =  u(-l)  =  0. 
Using  the  decomposition  method  [1],  let  L  a  d^/dx^  and  Lu  =  g  +  kxpu. 
Operating  with  L~\  we  have  L”\u  =  L~^g  +  L  ^kxpu.  Then 

u  =  C1  +  C2X  +  g*2/2  +  L 

00  o  -Id 

Let  u  *  l  un  with  uQ  =  c-j  +  cgx  +  gx  /2.  Then  um+1  *  L  kxHum  with 

n=0 

m  ^  0.  Thus  u  =  l  (L"^kxp)muQ  or  u  =  l  (l'^kxp)mc-j  +  l  (L  \xp)mc2x  + 
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l  (L  ^kxp)mgx2/2,  and  finally  u  =  c^<t^(x)  +  c2<*>2^x^  +  p(x)  where 


Mx)  =  l  kmxmp+2m/(mp  +  2m  -  1 ) (mp  +  2m) 
m=0 


$2(x)  =  l  kmxmp+2m+1/(mp  +  2m) (mp  +  2m  -  1) 
m=0 


r(x)  =  l  ( 1 /2 ) gkmxmp+2m+2/ (mp  +  2m  +  l)(mp  +  2m  +  2) 
m=0 


Since  u(l)  =  u(-l)  =  0  we  have  c^O)  +  c2«>2(l)  +  r(l)  =  0  and 
c-jt-j (-1 )  +  c2$2(-l)  +  r(-l)  =  0  or 


<t>](l)  4>20) 

^(-D  d>2  ( “1 ) 


-r(l) 


-r(-i) 


or  <|>C  =  r  or  C  =  <j>_1r  with 


M-l) 


-*2(1) 


^(DJ 

4>1(1)02(-1)  -  4>2(1  )4>-| (-1 ) 


C}  =  u2(i)r(-i)  -  <t>2(-i)r(i)]/Cd>1(i)4>2(-i)  -  ^2(1)^(-1)] 


C2  =  C4», (-1  )r(i )  -  <i>1(l)r(-i)]/[4>1(i )<|>2(-1 )  -  ♦2(i)*1(-i)] 


and  the  complete  solution  can  now  be  determined. 


Suppose  in  the  above  example,  we  let  k  =  40,  p  =  1,  g  =  2.  Thus 


we  consider  the  equation 


d2u/dx2  -  40xu  =  2 
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with  u ( — 1 )  =  u(l)  =  0.  This  is  the  one-dimensional  case  of  the  elliptic 
2 

equation  A  u  =  f(x,y,z)  +  k(x,y,z)u  arising  in  problems  of  physics  and 
engineering. 

2  2 

Here  let  L  =  d  /dx  and  we  have  Lx  =  2  +  40xu. 

This  Is  a  relatively  stiff  case  because  of  the  large  coefficient  of  u, 
and  the  nonzero  forcing  function  yields  an  additional  Airy-like  function. 
Operating  with  L"^  yields 

u  =  A  +  Bx  +  L~^(2)  +  L~^(40xu). 


un  =  A  +  Bx  +  L_1(2)  =  A  +  Bx  +  x2 


and  let  u  *  £  u  with  the  components  to  be  determined  so  that  the  sum  is  u. 
n=0  n 

We  identify  un+^  =  L"^(40xun).  Then  all  components  can  be  determined,  e.g 
u1  =  (20/3)Ax3  +  ( T  0/3 ) Bx4  +  2x5 


u2  =  (80/9)Ax6  +  (200/63)Bx7  +  (10/7)x8, 


etc.  An  n-term  approximant  d>  =  T  u.  with  n  =  12  for  values  of  x  is 

"  i=0  1 

given  in  the  table  below. 


n-term  approximant  <f>n  for  n 


-0.135649 

-0.113969 

-0.083321 

-0.050944 
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These  easily-obtained  results  for  a  12  term  approximation  are  exact 


to  seven  digit  accuracy. 

The  solution  Is  found  just  as  easily  for  nonlinear  versions  without 


linearization  by  simply  using  the  An  polynomials  for  the  nonlinear  terms. 
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APPROXIMATION  OF  NONLINEARITIES  USING  Ap  POLYNOMIALS 

Using  the  An  polynomials,  there  is  no  need  for  mathematically  inadequate 
and  physically  unrealistic  approximations  or  linearizations.  Thus,  if  the 
modeling  retains  the  inherent  nonlinearities,  we  may  expect  solutions  con¬ 
forming  much  more  closely  to  actual  behavior. 

Consider  simple  nonlinear  operators  not  involving  differentials,  i.e., 
of  the  form  Ny  =  f(y). 

For  polynomial  nonlinearities  the  An  being  sums  of  various  products 
of  the  y^  up  to  i  =  n  can  be  written  in  symmetrized  form.  Thus,  if  we 

oo 

have  Ny  =  y2  =  J  An,  AQ  =  yg,  A]  ■  2y()y1,  A?  =  y2  +  etc.;  but 

we  can  write  this  as  AQ  =  y^,  A]  =  yQy1  +  y^,  A2  =  y^  +  y]y]  +  y^Q, 
etc.,  i.e.,  the  first  subscript  goes  from  0  to  n,  and  the  second  is 
chosen  such  that  the  sum  of  subscripts  is  n. 

In  extending  to  any  analytic  function  f(y)  we  define  hn(y)  =  dnf/dyn 
and  write  a  convenient  heuristic  rule  which  users  may  find  more  convenient 
than  mathematical  derivations  - 

An  =  J,  c(v.n)hv(y0) 

For  example, 

A3  =  c(l,3)h1  +  c(2,3)h2  +  c(3,3)h3 

The  c(v,n)  can  be  obtained  by  simply  asking  how  ma-y  combinations  of  v 
integers  add  to  n.  Thus  c(v,n)  will  mean  the  sum  (from  1  to  v)  of  the 
products  of  v  of  the  y^  terms  whose  subscripts  add  to  n.  To  get  c(2,3), 
we  see  two  integers  can  add  to  3  only  if  one  integer  is  1  and  the  other  is 
2  (if  zero  is  excluded).  Hence,  we  write  c(2,3)  =  y^.  To  get  c(l,3). 


the  coefficient  of  h^(y0),  we  have  one  and  its  subscript  must  be  3, 
hence  c(l,3)  =  y3.  What  about  c(3,3),  the  coefficient  of  h3(yQ)?  Now  we 
need  3  factors  y^  with  subscripts  summing  to  3,  hence  each  subscript  must 
be  1  and  c(3,3)  =  y^yiy-j  *  y-j.  This  is  not  quite  right,  and  we  add  another 
heuristic  rule.  If  we  have  repetitions  of  subscripts  we  divide  by  the 
factorial  of  the  number  of  repetitions.  Then,  c(3,3)  =  (1/3! )y^ .  We  now  have 

A3  =  h1(y0)y3  +  h2^0^yly2  +  My0^1/3!)yl 

To  write  Ag,  for  example,  we  need  the  coefficients  for  the  terms  hv(yQ)  for 
v  from  1  to  6.  The  coefficient  of  hg  must  involve  six  integers  adding  to 
6  or  y®  hence  the  coefficient  of  hg(yg)  is  (l/6!)y^.  What  about  the 
coefficient  for  h^y^)  in  Ag  or  v  =  2,  n  =  6?  Clearly  we  need  two  integers 
that  sum  to  6.  These  are  (1,5),  (2,4),  and  (3,3).  Thus,  the  coefficient  c(2,6) 
is  (1/2! )y3  +  y2y4  +  y^g.  The  terms  involve  nV==1  yk  with  k.  =  n, 
and  if  we  have  j  repeated  subscripts,  we  divide  by  j!: 

A0  =  ho(yo} 

A1  =  hl(y0)yl 

A2  =  Wy2  +  h2^y0^1/2^yl 

A3  =  hl(y0)y3  +  h2(y0)yly2  +  h3<y0)(1/3’)yl 

A4  =  My(Py4  +  My0^1/2:)y2  +  yly3^ 

+  h3(y0)(l/2!)y2y2  +  h4(yQ) (1/4! )yA 

A5  *  hl(y0)y5  +  hZ{y0Ky^3  +  yly4] 

+  h3(y0)[yi(l/2!)y2  +  (l/2!)y2y3] 

+  h4(y0)(l/3!)y3y2  +  hg(y0) <1/5! )yf 
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>  l.t  li 


1  '*•  .**1 


it,  tl.  l>. 


A6  =  hl(y0)y6  +  h2(y0)[(1/2!)y3  +  y2y4  +  yly5] 

+  h3(y0)[(l/3:  )y|  +  y}yg3  +  (l/2!)y2y4] 

+  h4(y0)[(l/2:)y2(l/2!)y2  +  (l/3!)yfy3] 

+  h5(y0) (1/4! )y^y2  +  h6(y0)(i/6: )yf 

A7  =  hl(y0)y7  +  My0)Cy3y4  +  y2y5  +  yly6] 

+  h3(y0)[(l/2:)y2y3  +  y7(l/2:)y|  +  y^  +  (l/2:)yfy53 

+  h4(y0)Ly1(l/3!)y|  ♦  (1/2 :)y*y2y3  +  (l/3!)yfy43 
+  h5(y0)C(l/3!)y3(i/2:)y|  +  (l/4!)y}y3] 

+  Vyo)(1/5:>yiy2  +  h7(y0)(l/7:)y( 

A8  =  hl(y0)y8  +  h2^y0)t(1/2!jy4  +  ^5  +  y2y6  +  yly7] 

+  h3(y0)[y2(l/2!)y2  +  (1/2 !  )y2y4  +  y^  +  y^  +  (l/2:)y2y6] 
+  h4(y0)[(l/4: )y^  +  y](l/2:)y2y3  +  (1/2! )yf (1/2! )y2 
+  0/2:)y2y2y4  +  (1/3!  )yfy5] 

+  h5(y0)[(l/2:)y2(l/3!)y^  +  (l/SJJy^  +  (l/4!)y}y4] 

+  h6(y0)[(l/4:)yA(l/2:)y2  +  (l/5!)y*y3] 

+  h7(y0)(V6:)y*y2  +  hg(y0)(l/8!  )yf 
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Ag  =  Vy0)y9  +  h2(y0[y4y5  +  y3y6  +  y2y7  +  y^g] 

+  h3(y0)[(l/3!)y|  +  y2y3y4  +  (l/2!)y*y5  +  yi(l/2!)y* 

+  ylV5  +  yly2y6  +  (1/2!>yly73 
+  h4(y0)[(l/3:)y^y3  +  y}yz{V2\)y\  +  y7(l/2:)y|y4 

+  (1/21  )y7y3y4  +  (l/2I)yfy2y5  +  (l/3!)y^6] 

+  h5(y0)[yl(1/4*)y2  +  (1/2! )yf (1/2! )y*y3 
+  0/3!  )y^(l/2!  )y3  +  <1/3!  Jy^  +  (l/4!)y4y5] 

+  h6(y0)[O/3:  )yf  0/3!  )y|  +  (l/4!)y|y2y3  +  (l/5!)y*y43 
+  h7(y0)[(l/5!)yf(l/2!)y|  +  (1/6! )yfy33 
+  h8(y0)(l/7!)y(y2  +  h9(y0) (1/9! )y* 

*10  =  hl(y0)y10  +  h2(y0)C(1/2!)yl  +  y4y  6  +  y3y7  +  y2y8  +  yly9] 
+  h3(y0)C(l/2! )y|y4  +  y2(l/2!)y4  +  y2y3y5  +  (l/2!)y2y6 

+  Wy4y5  +  yly3y6  +  yly2y7  + 

+  h4(y0)C(l/2! )y|(l/2!)y2  +  (1/3! )y|y4  +  y1(l/3!)y^ 

+  yiy2y3y4 +  y^yz')y\ys  +  0/2!)yf(i/2:)y| 

+  (l/2!)yfy3y5  +  (1/2!  Jy^  +  (l/3!)yfy73 
+  h5(y0)[(l/5!)y|  +  y1(l/3!)y^y3  +  (1/2! )yfy2(l/2! )y\ 

+  (l/2!)yf(l/2!)y|y4  +  (1/3!^^ 
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+  O^Oy^y^  +  (l/4:)y^yg] 

+  h6(y0)C(l/2! )yf (1/4! )y|  +  (1/31 )yf (1/2! )y|y3 
+  (l/4!)y^(l/2:)y|  +  (l/4l)y*y2y4  +  (l/5!)y^y5] 

+  h7(y0)C(i/4l  )y^d/3!  )y|  +  (l/sOy^  +  d/6!)yfy4] 
+  hs(y0)C(l/6! )yf (1/2! )y|  +  (1/7!  Jy^,] 

+  h9(y0)(l/8!)yfy2  +  h10(yQ) ( 1/70! )y]° 


Recent  work  has  established  that  £  A  forms  a  generalized  Taylor  series 

n=0 

about  the  function  uQ(x). 


Computational  Time:  In  numerical  solutions  of  physical  problems,  it  is 
common  to  make  computations  at  discrete  space  or  time  intervals.  Computer 
methods  are  based  on  changing  continuous  problems  to  discrete  problems. 

Thus,  in  solving  a  differential  equation,  one  must  solve  the  equation  at 
each  point  of  time. 

Since  these  points  must  be  close  together  to  approximate  the  total  solu¬ 
tion,  massive  computations  are  needed  and  the  resulting  numerical  printouts 
yield  little  insight  into  dependences.  Another  drawback  of  these  methods  is 
the  fact  that  they  involve  linearized  approximations  in  each  interval.  As  the 
mesh  is  made  finer  to  increase  accuracy  the  number  of  computations  can  become 
enormous. 

For  example,  solving  even  the  simple  linear  equation 
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with  y(0)  s  1  and  y * (0)  =  0  requires  a  thousand  solutions  of  the 
above  equation  or  8,000  actual  computations  to  get  the  value  of  y 
every  .01  seconds  for  10  seconds.  For  one  minute  this  means  6,000  solu¬ 
tions  or  48,000  actual  computations.  Note  that  if  one  considers  an  equation 
in  two  independent  variables  x,t  we  require  small  steps  in  both  variables. 
Then,  the  computations  can  go  up  by  several  orders  of  magnitude. 

One  can  easily  visualize  a  problem  leading  to  a  billion  coupled  difference 
equations  to  solve.  For  a  simple  scalar  elliptic  equation,  we  have  one 
unknown  at  each  mesh  point.  For  more  complex  problems,  there  can  be  many 
unknowns  at  each  mesh  point  and  the  resulting  systems  of  difference  equations 
(instead  of  being  linear  as  in  the  previous  case)  may  be  nonlinear, 
time-dependent,  and  very  large  (inclusion  of  stochastic  coefficients,  etc., 
is  still  another  matter).  To  solve  massive  systems,  iterative  procedures 
are  used  to  solve  simpler  systems,  then  substitution  to  get  "residuals" 
and  repetitions  of  the  process  to  produce  corrections  and  until  the  error 
is  (or  i s  felt  to  be)  within  tolerable  limits.  To  get  accuracy  the  mesh 
must  ome  very  fine  and  computations  required  finally  exceed  any 
conceivable  computer  capability  for  complicated  equations  in  x,y,z,t. 

Solution  by  the  decomposition  method,  on  the  other  hand,  is 
continuous ,  analytical,  and  requires  no  discretization;  it  corresponds 
to  the  results  obtained  by  an  infinite  number  of  computations  by 
discretization  -  no  linearization  is  involved ,  and  the  solution  is 
accurate.  If  variable  coefficients,  several  independent  variables, 
and  nonlinearities  are  involved,  the  decomposition  method  is  clearly 
preferable.  In  the  case  of  stochastic  equations,  the  decomposition 
method  is  particularly  appropriate;  it  requires  no  perturbative  or 
truncation  methods  or  apriori  assumptions  of  special  behavior.  Computer 


a?  *»•  ■ 


results,  on  the  other  hand,  are  not  correct  when  stochastic  processes 
are  discretized. 

Supercomputers  are  developing  rapidly  because  of  urgent  need  in 
meteorology,  fluid  dynamics,  fusion  research,  intelligent  missile  guidance, 
and  weapons  design.  In  fluid  dynamics,  for  example,  they  are  considered 
essential  for  the  solution  of  the  equations  which  are  relevant  to 
turbulence,  internal  waves  in  the  ocean,  and  future  development  of 
hypersonic  flight  vehicles  and  engines.  Supercomputers  are  also  essential 
for  VLSI  devices,  seismology,  reservoir  modeling,  bioengineering,  and 
studies  of  the  national  economy.  The  decomposition  method  now  offers  the 
possibility  of  substantial  decrease  in  computational  time  over  current 
supercomputers. 

To  solve,  for  example,  the  relevant  problem  of  a  space  shuttle  for 
eventual  single-stage  flight  to  a  space  station  or  servicing  of  space 
structures,  a  three-dimensional  mesh  is  generated  which  discretizes  the 
system  of  nonlinear  partial  differential  equations  into  a  million,  a 
hundred  million,  or  perhaps  a  billion  coupled  difference  equations  in  as 
many  unknowns.  One  begins  to  see  then  the  tremendous  data  handling 
problem,  the  necessity  for  improved  algorithms  and  the  need  for  still 
greater  computational  speed.  We  may  also  have  many  unknowns  at  each  point, 
and,  as  we  have  pointed  out,  the  system  nonlinearities  and  random 
fluctuations  need  to  be  taken  into  consideration.  Since  usually  solutions 
are  iterative  -  first  solving  an  approximation  to  the  original  system  of 
differential  equations  and  then  improving  the  solution  by  repeated 
substitution  of  each  new  solution  -  parallel  processing  is  complicated  by 
the  difficulty  of  partitioning  the  work  so  each  processor  can  work 
independently.  This  is  being  pursued  by  many  ingenious  ideas  necessitated 
by  the  method  of  discretization. 
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In  all  such  problems  we  need  to  be  able  to  solve  coupled  systems 
of  nonlinear  (and  generally  stochastic  as  well)  partial  differential 

★ 

equations  with  complex  boundary  conditions  and  possible  delayed  effects. 

These  systems  are  linearized  and  discretized  (and  the  stochastic  aspects 
either  ignored  or  improperly  dealt  with)  so  the  various  numerical 
approximation  methods  can  be  used.  This  requires  faster  and  faster  super¬ 
computers  to  do  these  computations  in  a  reasonable  time.  Fifth  generation 
computers  are  being  considered  to  operate  at  speeds  up  to  1000  megaflops 

Q 

(a  gigaflop)  or  10  operations  per  second  or  even  higher. 

Unfortunately  the  further  developments  in  supercomputers  can  quite 
possibly  still  give  wrong  answers  because  even  a  single  one-dimensional 
nonlinear  differential  equation  without  stochasticity  in  coefficients, 
inputs,  and  boundary  conditions  -  let  alone  vector  partial  differential 
equations  in  space  and  time  with  nonlinear  and/or  stochastic  parameters 
arising  in  control  theory  for  space  structures. 

A  supercomputer  is,  after  all,  a  fast  adding  machine,  and  its  computa¬ 
tional  accuracy  Is  dependent  on  the  sophistication  of  the  mathematical  methods 
programmed  into  it.  Typical  calculations  consider  millions  of  discrete 
time  intervals  made  small  enough  so  trajectories  between  them  can  be  taken 
as  low-order  polynomials,  e.g.,  quadratics.  If  stochasticity  is  involved, 
then  Monte-Carlo  methods  are  used,  which  insert  randomness  but  not  the 
properly  correlated  randomness  which  is  present  in  the  physical  problem. 

When  one  studies  airflow  about  aircraft  surfaces,  computations  are 
made  at  tens  of  millions  of  points,  and  it  is  felt  that  increasing  the 


Delayed  effects  are  not  discussed  here  but  are  dealt  with  in  [2j. 
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volume  of  computation  to  the  limit  in  an  ultimate  extrapolation,  super¬ 
computers  will  yield  complete  accuracy.  Not  only  does  this  Ignore 
stochasticity,  it  ignores  the  sensitivity  of  nonlinear  stochastic  systems 
to  very  slight  changes  in  the  model  -  in  fact,  to  changes  essentially 
undeterminable  by  measurement. 
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DECOMPOSITION  TREATMENT  OF  NONLINEARITY 
COMPARED  WITH  LINEARIZATION 


Nonlinear  equations  arise  in  every  area  of  application,  and  the  correct 
solution  of  dynamical  systems  modeled  by  nonlinear  ordinary  differential 
equations,  systems  of  differential  equations,  partial  differential  equations, 
and  systems  of  partial  differential  equations  is  vital  to  progress  in  many 
fields. 

In  order  to  make  these  equations  tractable,  it  is  quite  common  to 
linearize  equations  or  assume  "weak"  nonlinearity,  etc.,  because  adequate 
methods  simply  have  not  been  otherwise  available.  The  practice  of 
approximating  a  nonlinear  function  with  a  linearized  version  arose  from  the 
need  to  make  equations  tractable  by  simple  analysis,  since  numerical 
solutions  from  computers  can  have  drawbacks  in  that  functional  relationships 
are  difficult  to  see  in  numerical  printouts  and,  in  many  cases,  computation 
can  be  excessive,  and  methods  of  analytical  solution  of  nonlinear  equations 
have  been  generally  inadequate.  It  is  known,  of  course,  that  the 
linearized  solution  can  deviate  considerably  from  the  actual  solution  of 
the  nonlinear  problem  and  that  linearization  procedures  require  proof 
that  the  solution  is  valid.  For  example,  writing  x"  =  a  sin  x  in  the 
form  x"  =  ax  requires  a  priori  proof  that  x  is  sufficiently  small. 

Usage  of  linearization  has  become  rather  standardized;  however,  solution 
of  the  actual  nonlinear  form  is  clearly  preferable  to  a  linear  approximation. 

The  decomposition  method  has  substantially  improved  our  ability  to 
solve  a  wide  class  of  nonlinear  and/or  stochastic  equations.  It  is  now 
possible  to  obtain  very  accurate  and  verifiable  solutions  of  nonlinear, 
or  even  nonlinear  stochastic  equations  for  all  of  the  above  types  even  if 
nonlinear,  stochastic  or  coupled  boundary  conditions  are  involved. 


Until  now,  linearization  and  perturbation  have  been  essential  procedures 


and  one  is  reluctant  to  give  up  such  a  convenient  analytical  tool  as 
superposition.  However  even  with  the  increases  in  computer  speed,  the 
complexity  of  some  of  the  problems  demands  new  techniques  to  save  computa¬ 
tional  time.  This  will  be  an  important  advantage  of  decomposition. 

Exact  linearization*  is  possible  for  some  nonlinear  equations  so  that 
a  convenient  check  can  be  made  of  the  decomposition  solution.  As  an  example, 
consider  the  nonlinear  equation  for  <j>(x,t) 

(j)^  +  <j)^  +  a$  +  <j>*"  =  0 

with  specified  conditions.  The  transformation  $  =  l/<j;(x,t)  leads  to 
the  linear  equation 

i|>t  +  ~  =  1 

and  conditions  specified  on  <J/.  This  is  now  a  linear  equation  which 
presents  no  difficulty  if  the  conditions  on  are  specified.  However, 
the  nonlinear  equation  can  be  solved  directly  as  follows.  Using  decomposition 
we  write 

Lt <J>  +  Lx<|>  +  a<j>  +  <j>2  =  0 
00 

Let  N<j>  *  a<j>2  5  I  A  where 
n=0  n 

A0  =  *o 

A]  =  2Vl 

★ 

By  exact  linearization  we  refer  to  the  cases  where  it  is  possible  to 
use  transformations  of  dependent  and  independent  variables  to  transform  a 
nonlinear  equation  into  a  linear  equation. 
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A2  s  $1  +  2<*)o^2 

^3  =  +  2<t>0^3 

Solving  for  L.<j>  and  for  L  <$> 

t  X 

=  -Lx<#>  -  a<t>  -  <t>^ 

Lx<j>  =  “Lt*  '  a<J>  "  4>2 

Thus 

<p  =  4>(x,0)  -  L*1LX«#>  -  aL"1*  -  L’V 
<P  =  $(0,t)  -  L"1^  -  aL"1*  -  L^V 

Adding 

4>  =  (1/2){<|>(x,0)  +  <j>(0,t)  -  (L"\x  +  L"1Lt)4) 

~a(Lt1  +  Lx1),p  +  (Lt1  + t;V> 

We  define 

<j>0  =  (1/2){4>(x,0)  +  4>(0  9 1) } 

00  OO 

substitute  <f>  «  l  <t>  and  <p2  =  J  An  to  obtain 
n=0  n  n=0  n 


Vi '  -0/2>a-t\  -  L;'Lt)*„ 

-  0/2WL;1  +  L;1  }<frn  -  (1/2){L^1  ♦  L;1  }An 

for  n  >  0  so  that  all  components  can  be  determined.  The  solution,  of 
course,  will  be  identical.  However  most  equations  cannot  be  exactly 
linearized  by  transformations,  hence,  the  decomposition  technique  is  much 


more  useful. 
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In  the  analysis  of  systems  it  is  common  to  suppose  that  everything  is 
linear  or  sufficiently  close  to  linear  so  that  linearized  analyses  will 
be  adequate.  Thus  in  a  mechanical  system  assumed  to  be  linear, 
displacements  and  accelerations  are  proportional  to  forces.  Thus  we  write 
output  y(t)  as  proportional  to  input  x(t)  or  y  =  kx  where  k  is  a 
constant  independent  of  t  or  x.  Now  suppose  that  the  system  deviates 
only  slightly  from  linearity  by  adding  a  small  term  which  is  nonlinear,  e.g., 
y  =  kx  +  ex2.  If  we  now  consider  the  input  x  =  A  cos  w-jt  +  B  cos  u^t,  we 
get  not  only  the  linear  term 

K[A  cos  +  B  cos 

but  also 

k  [A2  cos2  w-jt  +  B2  cos2  w2t  +  2AB  cos  u^t  cos  ^t] 

The  first  two  of  these  produce  constant  terms  and  second  harmonic  terms 
as  before.  Also  sum  and  difference  frequencies  arise  from  the  cross 
product  term 

2AB  cos  o^t  cos  Ugt  =  AB[cos(w-j  +  u>2)t  +  cos(uvj  -  o)2)t] 

Thus,  a  nonlinear  system  produces  new  effects  not  present  in  linear 
systems:  these  effects  are  proportional  to  e  (and  to  products  of 
amplitudes  A2,  B2,  or  AB).  Clearly  then,  if  e  is  not  small,  such 
effects  become  important. 

Solutions  are  generally  carried  out  only  under  the  assumption  that 
e  is  small  so  that  perturbation  theory  will  be  applicable,  i.e.,  when 
we  consider  a  "slightly  nonlinear"  or  a  "weakly  nonlinear"  system. 

Since,  in  general,  solutions  of  nonlinear  equations  are  made  by 
linearizing  the  equations,  it  is  natural  to  ask  what  the  effect  of 


linearization  is  on  the  actual  solutions.  Let  us  consider  the  general 
nonlinear  example 


Ly  +  Ry  +  Ny  =  x(t) 


where  L  is  the  invertible  linear  operator,  R  is  the  remaining  linear 
operator,  and  Ny  is  the  nonlinear  term.  We  have 


Ly  =  x  -  Ry  -  Ny 


y  =  <j>  +  L_1x  -  L_1Ry  -  L‘]Ny 


where  L <p  =  0.  Assume  the  solution  is  given  by  y  =  [  yB  with 

n=0  n 

yg  =  <j>  +  L~^x  identified  as  the  first  component  of  the  sum.  Assume  also 

oo 

the  decomposition  of  the  nonlinear  term  Ny  into  £  An  where  the  An 

are  generated  for  the  specific  Ny.  Then  the  components  after  yQ  are 
determinable  in  terms  of  yQ  as 


*i =  -L'lRy0  -  L_1W 

y2  =  -L'1Ry1  -  L‘%  (yg^ ) 

Yg  =  ”L  Ry2  ~  L  A2(yg,y.|  ,y2) 


yn  =  "L  lRyn-l  "  L  lAn-l^y0”-,,yn-l^ 


or,  equivalently, 


=  ( -l”1 R)y0  -  l_1a0 

y2  =  (-L_1R)2y0  -  (-L'1R)L_,A0  -  L“1A1 

y3  -  (-L-1R)3y0  -  (-L'1R)2L'1A0  -  (-L_1R)L’1A1  -  L_1A2 
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yn  -  (-L_1R)ny0  -  l  (-L‘1R)n'1"v  -  L']A 

v=0 

00 

for  n  _>  1.  The  solution  is  y  =  l  y  or 

n=0  n 

y  =  l  {-L’\)\ 
n=0  u 

00  n-1  ,  ,  , 

-  I  l  (-L-1R)"-1-v  -  L-’A, 
n=l  v=0 

It  has  been  shown  by  Adomian  that  An(y0,y.| » •  •  *  »yn)  reduces  to  yn  if 
Ny  =  f(y)  =  y.  Then  the  solution  is 


y  =  I  (-L_1R)ny0 
n=0 

oo  n-1 

-  I  l  (-l-'r) 

n=l  v=0 


-lnNn-l-v_L~ly 


i.e.,  the  solution  corresponds  now  to  the  equation  Ly  +  Ry  +  y  =  x 
which  yields 

y  =  y0  -  L_1Ry  -  L-1y 

with 

yi  =  "L’lRyo  "  L"1yo 
y2  =  -L’V,  -  l‘\ 

=  (-L_1R)2y0  -  (-L~1R)L_1y0  -  l\Q 
y3  =  "L  lRy2  "  L  1yl 

=  (-L_1R)3y0  -  (-L“1R)2L~1y0  -  (-L-1R)L_1y0  -  L"1y1 
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The  result  on  the  solution  of  replacing  Ny  =  f(y)  by  y  can  be 
determined  by  plotting  the  nonlinear  and  the  linearized  results  for  a  specific 
f(y)  as  we  will  see  in  the  following  example.  Similarly  replacing  f(y) 
by  another  perhaps  more  sophisticated  linearization  is  seen  by  simply 
calculating  the  An  for  the  linearized  function  replacing  f(y).. 


Example:  Exponential  nonlinearity:  Let  us  consider  a  simple  nonlinear 
ordinary  differential  equation  with  an  exponential  nonlinearity 

dy/dx  +  ey  =  0 
y(o)  =  1 

In  our  usual  standard  form  (1983)  this  is  written 


Ly  +  Ny  =  x 


with  L  =  d/dx  and  Ny  =  ey.  We  solve  for  Ly,  i.e.,  Ly  =  -Ny  hen 
write  L  \y  =  -L  ^Ny  with  L~^  defined  as  the  integration  over  x. 
The  left  side  is  y  -  y (0)  hence 

y  =  y(0)  -  L_1Ny 


The  nonlinear  term  Ny  =  ey  is  replaced  by  Y  A„  where  the  A„ 

n=0  n 

can  be  written  as  ^(e^)  to  emphasize  that  they  are  generated  for 
this  specific  function.  Thus, 


y  =  y (0)  -  L'1  l  A  (ey) 
n=0  n 

Now  the  decomposition  of  the  solution  y  into 
term-by-term  identification 


00 


l  yr 

n=0  r 


leads  to  the 


* 


See  reference 
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y0  *  y(o)  =  i 


y,  -  -l'1ac 


y2  *  -L  'a, 


-L~ 'A. 


The  An(ey)  are  given  by 


v" 


A1  =  y^ 


A2  =  (yf/2  +  y2)e> 


A3  =  (y^/6  +  y]y2  +  y3)ei 


y-,  =  -ex 


y2  =  e2x2/2 


3„3  i-i  i 


y-,  =  -e  x  /3 


y  =  1  +  l  (-1 )nenxn/n! 
n=l 


which  can  also  be  written  y  -  1  -  ln\_  1  +  ex]  for  x  <  1 
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Suppose  we  approximate  the  e^  term  by  1  +  y,  dropping  all  terms 
of  the  series  for  e^  except  the  constant  and  the  linear  term.  The 
differential  equation  becomes 

dy/dx  =  -(1  +  y) 

Then 

y  =  y ( 0)  -  L"1 (1  +  y) 

=  1  -  L-1 [1 ]  -  L_1  l  y 

n=0  n 


so  that 


y0  =  l  -  x 


y.  »  -L"  (1  -  x) 


y„  =  -L~yn.,  ('ll) 


Since  y  is  the  sum  of  the  components,  we  have 


y  =  ?  {(-1 )nxn/n!  +  (-l)n+1xn+1/(n+l)!} 

f 

n=0 

y  =  e'x  +  (e"x  -  1)  =  2e"x  -  1 

§ 

We  will  identify  this  linearized  solution  as  y  and 

compare  with  the 

w*. 

solution  y  of  the  nonlinear  equation.  The  results 

are  given  in  the 

v» 

table 
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A  COMPARISON  OF  ACTUAL  AND  LINEARIZED  SOLUTIONS 
(EXPONENTIAL  NONLINEARITY) 

Linearized  Solution 

Actual  Solution  y(x)  y0(x) 


1.0000 

.7595 

.5658 

.4036 


.1417 


V*> 

1.0000 

.8097 

.6375 


.3406 


Example  (Anharmonic  oscillator):  As  an  example  of  a  vibration  problem, 

2  2  2 

consider  now  the  well-known  anharmonic  oscillator  equation  d  0/dt  +  k  sin  0 
for  0(0)  =  y  =  constant  and  0'(O)  =  0.  Using  the  decomposition  method, 
the  solution  is  found  to  be 

0(t)  =  y  -  [(kt)2/2!]sin  y  +  [(kt)4/4!]sin  y  cos  y 


-  [(kt)6/6!](sin  y  cos2  y  -  3  sin3  y]  + 


which  becomes 


6(t)  =  Y[1  -  (k2t2/2! )  +  (kt)4/4:  -  •••] 

in  the  linearized  case,  i.e.,  for  "small  amplitude"  motion  which  offers 
interesting  comparison  in  order  to  determine  when  the  smallness  assumption 
is  inappropriate  because  of  large  amplitude  vibrations. 
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Example:  Hyperbolic  sine  nonlinearity:  Consider  the  equation 
du/dt  =  k  sinh  u/a 
where  u{0)  =  c  >  0  for  t  >  0. 

If  we  assume  that  we  can  approximate  sinh  u/a  «  u/a,  we  have 

oo 

du/dt  -  ku/a  =  0.  Now  solving  by  decomposition  with  u  =  \  u  , 

n=0  n 

L  =  d/dt,  and  L-^  as  the  definite  integral  from  0  to  t 


Lu  -  ku/a  =  0 


L_1Lu  =  L_1 (k/a)u 


u  =  u(0)  +  L~  (k/a)u  =  l  u 

n=0  n 

Uq  =  u(0)  =  c 


u-j  =  (k/a)ct  =  ckt/a 
u2  =  c(kt/a)2/2: 
u3  =  c(kt/a)3/3! 


um  = 


i.e.,  u  =  ce 


To  solve  the  original  equation  with  sinh  u/a 


u  =  l  u  -  c  +  kL"1  l  A 
n=0  n  n=0  n 

where  the  An  are  generated  for  Nu  =  sinh  u/a.  These  are  given  by: 
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is 
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s1nh(uQ/a) 

(Uj/a)  cosh(Ug/a) 

(u2/a)  cosh{Ug/a)  +  (l/2!)(iiZ/a2)  sinh(uQ/a) 
(u3/a)  cosh(Ug/a)  +  (u^/a)(u2/a)  sinh(uQ/a) 

+  (1/3! )(u^/a3)  cosh(Ug/a) 

(u4/a)  cosh(uQ/a)  +  [(1/2! ) (u|/a2) 

+  (u1/a)(u2/a)]  sinh(Ug/a) 

+  (l/2!)(u^/aZ)(u2/a)  cosh(Ug/a) 

+  (1/4! ) (u^a4)  s1nh(u0/a) 

(u5/a)  cosh(uQ/a)  +  [(u2/a) (u3/a) 

+  (Uj/aMu^/a)]  sinh(uQ/a) 

+  [(u1/a)(l/2! ) (u2/aZ) 

+  (l/2!)(uZ/a2)(u3/a)]  cosh(u0/a) 

+  (1/3! )(u^/a3)u2/a)  sinh(u0/a) 

+  (1/51 ) (u^/ot5)  cosh(uQ/a) 

(uga)  cosh(uQ/a)  +  [(1/2! )(u|/a2)  +  (u2/a)u4/a) 
+  (u-j/a)(ug/a)]  sinh(Ug/a) 

+  [(1/3! )(u~/a3)  +  (u,/ct)(u,/a)(u,/a) 


+  (1/2!)(uj/oT)(u4/a)]  cosh(Ug/a) 

+  [(1/2! )(u2/a2)(l/2! )(u2/a2) 

+  (l/3!)(u3/a3)(u3/a)]  sinh(u0/a) 

+  (l/4!)(u|/a4)(u2/a)  cosh(u0/a) 

+  (1/6! )(u®/a6)  sinh(uQ/a) 

c 

kL_1A0  =  kL"^[sinh  uQ/a]  =  kt  sinh  c/a 
kL-1  A-j  =  kL"1[(u1/a)  cosh(uQ/a)] 
kL"^[(kt/a)  sinh(c/a)  cosh(c/a)] 

(k2t2/2!a)  sinh(c/a)  cosh(c/a) 

kL  \  =  kL_1[(u2/a)  cosh(u0/a)  +  (l/2)(u2/a2)  sinh(uQ/a)] 
kL  ^[(k2t2/2)(l/a)2  sinh(c/a)  cosh2(c/a) 

+  (l/2)(k2t2)(l/a)2  sinh3(c/a)] 

(k3t3/3! )(l/a)2  sinh(c/a)[sinh2(c/a)  +  cosh2(c/a)] 

(k4t^/4! )(l/a)3[sinh(c/a)  cosh(c/a)][5  sinh2(c/a)  +  cosh2(c/a)] 
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thus,  the  correct  solution  is  the  sum  of  the  u  above  while  the  linearized 

n 

solution  is 

u  =  c[l  +  kt/a  +  (kt/a)2/2!  +  (kt/a)3/3!  +  •••] 

If  we  assume  a  *  c  =  k  =  1,  we  have  the  results  u  =  e*  in  the  linear 
case  which  can  be  compared  with  the  nonlinear  solution. 


ERROR  IN  LINEARIZED  SOLUTION 
(FOR  A  HYPERBOLIC  SINE  NONLINEARITY) 


t 

Linearized  Solution 

Nonlinear  Solution 

%  Error 

0 

1 

1 

0 

.1 

1.105170918 

1.127402502 

1.97% 

.2 

1.221402758 

1.278624737 

4.48% 

.3 

1.349858808 

1.462380247 

7.69% 

.4 

1.491824698 

1.693614375 

11.92% 

.5 

1.648721271 

2.001468676 

17.62% 

.6 

1.8221188 

2.456234584 

25.82% 

.7 

2.013752707 

3.325545159 

39.45% 

.75 

2.117000017 

4.512775469 

53.09% 

.76 

2.13827622 

5.121285511 

58.25% 

.77 

2.159766254 

6.939848656 

68.88% 

.7719 

2.163873711 

10.9022661 

80.15% 

.771936 

2.163951611 

14.69149181 

85.27% 

.77193683 

2.163953407 

20.34929139 

89.37% 

.7719368329 

2.163953414 

28.3241683 

92.36% 

.7719368330 

2.163953414 

00 

100.  % 

So  we  see  that  the  error  due  to  linearization  approaches  100%  in  this 
simple  example.  Linearization  of  nonlinearities  in  modelling  of  rapidly 
maneuvering  space  structures  could  conceivably  result  in  similar  inaccuracies 
and  should  be  avoided. 
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and  since  AQ  =  y"1 ,  A1  =  -y^y-, ,  A2  =  -yg2y2  +  ygM’ 

yi  -  f  yn]dt  =  f  (k  +  t2/2)_1dt 
1  Jo  u  0 

y,  =  (2/k)1/2tan_1[t/(2k)1/2] 


Let  us  consider  a  two-term  approximation  4>2  =  yo  +  yl 

00 

solution,  of  course,  is  J  yn.)  Then 

<j>9  =  k  +  t2/2  +  (2/k)1/2tan-1[t/(2k)1/2] 


(The  complete 
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Series  Solution  Convergence  Compared  with  Numerical  Integration:  The 
decomposition  method  converges  very  rapidly  in  most  cases:  <j>n,  the  n-term 
approximation,  is  accurate  for  quite  low  values  of  n.  To  emphasize  this 
point  we  now  consider  only  a  two-term  approximation  in  the  following  example 
with  the  added  comment  that  additional  terms  are  extremely  easy  to  obtain. 

Consider  the  equation  dy/dt  =  t  +  y~^.  We  now  have  L  =  d/dt, 

Ny  =  -y'1,  x(t)  *  t.  Assume  y (0)  -  k  is  an  integer.  Then 

L_1Ly  -  L_1x  -  L_1Ny 

oo 

y  -  y(0)  -  L_1t  +  L_1  J  An 

n=0  n 

where  the  An  =  An(y_1). 

y  -  y(o)  +  L'1  l  A 
n=0  n 


y0  -  k  +  t  /2 
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The  following  table  compares  this  approximation  with  results  of  a  numerical 
integration  using  k  =  4.  With  only  a  little  more  effort  we  could  go  to  a 
higher  n  in  4>n  for  a  better  approximation,  which  is  necessary  since  the 
percentage  error  is  already  extremely  small.  The  worst  case  is  less  than 
0.4*.  However,  if  we  go  to  <j>j  we  find  the  worst  case  has  an  error  less 
than  0.02*  -  this  for  only  a  three  term  approximation.'  Thus,  we  have  very 
rapid  convergence. 
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COMPARISON  OF  DECOMPOSITION  AND  NUMERICAL  INTEGRATION 


Decomposition 

method 

4*1 


Numerical 

integration 
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HEAT  AND  DIFFUSION 


Typical  problems  of  this  type  are  initial-boundary  value  problems  over 
a  specified  region  and  with  a  given  initial  temperature.  Let's  look  at 
some  common  equations  involving  nonlinear  terms  which  may  arise  in  a 
particular  situation  -  all  of  which  are  easily  dealt  with. 

Burger's  equation  u^.  +  uux  =  uxx  represents  a  balance  between 
diffusion  and  convection.  If  the  uxx  term  vanishes,  we  have  a  simple 
diffusion  equation  -  a  scalar  hyperbolic  equation  u  =  u  .  Some  specific 

L  A  A 

conditions  must  be  specified  for  completeness,  e.g., 

u(x,0)  =  f(x)  for  -oo  <  x  <  0° 

u(0,t)  =  0 
and 

ux(0,t)  =  h(t) 

Another  possibility  is 

ut  +  uux  +  Xu  =  auxx  X  >  ° 

or  the  inhomogeneous  Burger's  equation 

"t 4  “« = 

for  0  <_  x  <_  i  and  0  <_  t  <_  T  with  u(0,t)  =  u(L,t)  =  0  and  u(x,0)  =  a(x) 
where  a(0)  =  a (L)  =0  and  a  >  0. 

Fisher's  equation  u^  =  uxx  +  u(l  -  u)  shows  still  another  nonlinearity 
u(l-u). 

For  a  general  nonlinear  form,  we  can  consider 
(3/3x)[k(u)3u/3x]  =  cp3u/3t 


39 


.v**'  a*  JtjP 


t  *  |  |  -*  | •  *  4’  1 


with  cp  a  constant  and  k(u)  a  specified  function.  If  k  =  1,  this 
becomes  the  simple  heat  equation  uxx  =  ut<  When  temperature  variations 
increase,  one  must  consider  k  dependent  on  temperature  u,  that  is, 
k  =  k(u)  so  that  the  equation  becomes  nonlinear.  Suppose,  for  example, 
k  =  kQ  +  k^u,  etc.  Finally,  we  can  have  equations  such  as 
ut  =  uxx  +  3f(u»ux)/3t  or  =  a2uxx  +  f(x,t),  etc.  All  of  the  foregoing 
are  easily  dealt  with  by  decomposition  with  the  An  polynomials  generated 
for  the  particular  nonlinear  terms. 


SOME  SPECIAL  CASES  DEPENDENT  ON  STRUCTURE  AND  MATERIALS 

1)  One-dimensional  heat  equation:  Consider  a  single  space  dimension  with  a 
specified  temperature  distribution  u(x,t)  described  by 


(3/3x)(k  3u/3x)  +  F(x,t)  =  cp(3u/3t) 


where  F  represents  a  heat  source  at  x  at  time  t.  In  the  mathematically 
simplest  case  where  k,  c,  p  are  assumed  to  be  constants,  we  can  write 
this  in  the  form 


iS 

m 
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ut  =  a2uxx  +  f(x,t) 


where  a d  =  k/cp  and  f(x,t)  =  (l/cp)F(x,t)  or  if  the  source  vanishes. 


simply 


ut  = 


often  referred  to  as  the  heat  conduction  equation.  We  can  consider,  for 
example,  a  homogeneous  rod  of  length  l,  thermally  Insulated  on  lateral 
surfaces  and  sufficiently  thin  so  that  at  any  time  t,  the  temperature 
u(x,t)  is  the  same  at  all  points  of  the  cross-section  at  x. 


•ViVt 


If  the  rod  is  inhomogeneous  so  that  k  =  k(x)  and  heat  exchange  takes 
place  with  the  surrounding  medium,  the  situation  is  described  by  an  equation 
of  the  form 

ut  -  a  uxx  +  au  =  f(x.t) 

with  a  =  h/cp,  h  a  heat-exchange  coefficient,  and  f(x,t)  = 
aT(x,t)  +  g(x,t)/cp  with  T  the  temperature  of  the  medium  and  g(x,t) 
the  density  of  the  heat  sources. 

2)  Three-dimensional  heat  equation:  In  the  case  of  three  space  dimensions 
and  time,  the  heat  flow  is  described  by  u(x,y,z,t)  and  the  equation 

cput  «  div(k  grad  u)  +  F 

where  k,  c,  p  are  functions  of  x,  y,  z.  More  generally  we  may  have  k 
dependent  on  u  and  space  coordinates,  e.g.,  in  the  case  of  large 
temperature  fluctuation.  This  is  the  more  physically  realistic  case  which 
we  can  apply  to  real  trusses  and  completely  amenable  to  solution  by  the 
decomposition  method. 

If  the  material  is  homogeneous,  we  can  write 
Ut  ’  a2(\x  +  uyy  +  UZZ>  +  F/c» 

2 

where  a  =  k/cp,  or, 

ut  *  a2V2u  +  f  f  =  F/cp. 

The  diffusion  equation  is  analogous  to  the  heat  conduction  equation. 
For  one  space  dimension  and  time,  we  can  write 


3/3x(D  3u/3x)  =  c(3u/3t) 


or  (Du  )  =  cuf.  If  the  diffusion  coefficient  D  is  constant,  we  have 

u.  *  a2u 
t  xx 

2 

in  the  same  form  as  before  with  a  =  D/c. 


3)  Boundary  conditions  for  Heat  -problems :  Differential  equations  do  not 
completely  specify  solution  uniquely.  We  must  also  prescribe  conditions 
on  the  solution.  Thus  to  complete  the  specification  of  the  problem,  we 
require  the  auxiliary  initial  and/or  boundary  conditions  as  well.  Thus  in 
the  heat  equation,  we  can  prescribe  the  temperature  at  the  ends  x  =  0 
and  x  =  l.  We  can  specify  u(0,t)  =  g(t)  where  g(t)  is  known  in  the 
time  interval  of  interest.  Or,  we  can  specify  the  heat  flow  at  one  end, 
e.g.,  3u(0,t)/3x  =  h(t).  Or,  we  can  have  more  complicated  conditions  such 

as 

3u(£,t)/3x  +  Au(£,t)  =  T(t) 

These  restrictions  apply  to  all  partial  differential  equations.  Such 
equations  have  many  solutions.  In  order  to  identify  a  particular  solution 
from  the  entire  possible  set  of  solutions,  we  need  supplementary  conditions 
similar  to  specification  of  initial  conditions  in  an  initial-value  problem 
described  by  an  ordinary  differential  equation. 

The  usual  situation  is  that  we  will  have  some  boundary  conditions 
specified  on  the  boundaries  of  the  space  described  by  the  coordinates 
x-j,  Xg»  x3  or  x,  y,  z  for  which  the  partial  differential  equation  is 
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applicable  and,  also.  Initial  conditions  at  a  particular  instant  t  =  0  or 
t  =  t0. 

If  we  know  the  temperature  distribution  at  t  =  0  and  the  temperature 
at  the  boundaries,  the  temperature  distribution  will  be  precisely  determined 
for  t  >  0.  If  the  time  interval  of  interest  is  very  long,  the  influence 
of  the  initial  condition  vanishes  and  the  distribution  will  be  determined  by 
the  boundary  conditions.  For  0  <_  x  <_  1,  we  might  specify  u(0,t)  =  f(t) 
and  u(£,t)  =  g(t).  For  0  <  x  <  ~,  we  need  only  u(0,t)  =  f(t). 

If  the  region  is  unbounded,  then  a  bounded  solution  is  uniquely 
determined  by  the  initial  conditions  alone,  i.e.,  by  u(0,u^ •  Thus 
for  the  bounded  case  we  might  have  conditions  such  as  u  at  t  =  0  and 
u  at  x  *  0  and  x  =  a. 

However,  if  for  very  large  a,  we  retain  the  conditions  at  t  =  0  and 

x  =  0,  but  now  u ( a )  approaches  zero  as  a  becomes  very  large. 

Use  of  the  given  conditions  is  simple.  Thus  in  a  second-order 

2  2 

differential  equation  with  no  forcing  term,  i.e.,  one  involving  d  u/dx  , 

00 

our  Uq  term  in  the  decomposition  u  *  £  up  is  given  by  u0  =  Y1  +  Y2X' 

n=0 

If  the  temperature  u  is  given  at  the  ends  x  =  0  and  x  =  a  of  the  bar, 
we  can  evaluate  y-j  and  in  terms  of  the  given  temperatures. 


4)  Solution  of  a  Two-dimensional  heat-flow  equation:  Consider  the  parabolic 

equation  utfV  +  u  =  kut  where  k  =  k(x,y).  Write  this  as 
xx  yy  t 

V  +  V  =  k(x>y)Ltu 

We  write  three  equations 

Lxxu  =  kLtu  -  Lyyu 

Lyyu  '  kLtu  ‘  Lxxu 

Ltu  =  k  tLxx  +  Lyy-^u 

Operate  on  the  first  with  L~],  the  second  L“J,  and  the  third  with  L,1 

xx  yy  i* 


to  get 


u  *  +  V  +  L«kLtu  -  LxxLyyu 

u  =  u(x,0,t)  ♦  Yjy  +  LyJkLtu  "  LyyLxxu 
u  =  u(x,y,0)  +  L^V1[Lxx  +  Lyy]u 
Adding  and  dividing  by  three 


u  =  (l/3){u(0,y,t)  +  u(x,0,t)  +  u(x,y,0)} 


+  (l/3){Ylx  +  Y2y} 

+  (1/3){[L~X  +  L'y]kLtu 
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Now  u  Is  replaced  by  l  u_  with  un  identified  as 

n-0  n  u 

uQ  =  (l/3){u(0,y,t)  +  u(x,0,t) 

+  u(x,y,0)  +  y-j*  +  Y2y) 


The  boundary  conditions  at  x  =  a  and  y  =  b  determine  y-j  »  Y2-  These 

vanish  if  a  and  b  recede  to  infinity.  Each  u„  after  un  is  now 

n  u 

determined  by: 


u  = 


0/3){[L^  ♦  L;}]kLtVl 

^•LxxLyy  +  LyyLxx-^un-l 
+  Lt  k  tLxx  +  Lyy3Un_1 


L  involves  two  differentiations  with  respect  to  x.  L  involves  two 
**  yy 

differentiations  with  respect  to  y.  l_t  involves  a  single  differentiation 
with  respect  to  t.  The  inverses  are  the  corresponding  integrations. 

oo 

Thus,  all  the  components  of  u  =  ][  u  are  determinable  and  again  we 

n=0  n 

see  that  the  decomposition  method  solves  such  equations. 


5)  Solution  of  a  two-dimensional  heat- flow  equation  with  a  heat  source  g: 
We  will  consider  a  bounded  rectangular  plate  with  temperature  specified  on 
the  edges  and  at  an  initial  time  t  =  0.  The  equation  to  be  considered  is 

uxx  +  “yy  -  k<*oOut  -  g 
or  equivalently 

92u/9x2  +  92u/9y2  -  k(x,y)  9u/9t  =  g 


45 


for  0  <  t  <  »  and  x,y  e  fi,  a  bounded  rectangular  plate  with  boundary 
T  with  temperature  distribution  specified  on  the  edge  x  e  [0,a],  and 
on  y  e  [0,b],  k  =  k(x,y),  given  u(x,y,0)|^  =  0  and  u(x,y,t)  =  f(x,y,t) 
for  x,y  e  r. 


Suppose  now  we  let  the  conductivity  k  be  constant  and  the  temperature 
T  be  constant  on  the  given  surfaces.  Then 


un  =  T  -  (Tx/a  +  Ty/b)/3 


If  a,b  -*•  then  uQ  =  T,  i.e.,  uQ  approaches  T  for  large  a,b 

U1  “  1/3fLxx  +  V“o  "  1//3^LxxLyy  +  LyyLxx-^u0 
+  1/3  L;V(:Lxx  *  Lyy]u0  -  (T/6)(x2  +  y2) 

oo 

We  know  Ug,  the  first  component  of  u  =  £  un  and  we  have  now  found  u-j 
in  terms  of  Ug.  All  other  components  are  found  in  the  same  way.  If 
we  replace  u-j  by  un+^  and  Ug  by  up  and  let  n  =  0,  we  have  the 
above  formula.  If  we  let  n  =  1,  we  have  in  terms  of  u-j.  Then 
let  n  =  2,3,...  to  get  others. 
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6)  Three-dimensional  general  case;  Let  us  consider  the  general  form 
where  k  may  be  a  constant,  space  or  time  dependent,  or  a  function  of  u. 
In  the  last  case,  we  have  a  nonlinear  equation.  The  nonlinear  terms  must 
be  replaced  by  the  sum  of  the  appropriate  Ap  as  specified  by  the 
decomposition  method 
3 

Y  (9/9x.)(k  9u/3x . )  -  cp3u/3t  =  g 
i=l  1  1 

where  k  is  the  coefficient  of  thermal  conductivity,  p  is  the  density, 
and  c  is  the  specific  heat  at  x^,  x^,  x^.  We  will  use  x,  y,  z  to 
avoid  confusion  with  the  subscripts  for  the  decomposition.  Let  Lx  =  3/3x, 
L =  9/9y,  Lz  =  9/3z,  Lt  =  9/9t.  Now  we  have 

LxkLxu  +  LykLyU  +  LzkLzu  -  cpLtu  =  g 


Solving  for  each  term  in  turn. 


LxkLxu  =  9  +  CpLyu  '  LykLyu  "  LzkLzU 


LykLyU  =  g  +  cpLtu  -  LxkLxu  -  LzkLzu 
LzkLzu  =  9  +  cpLtu  "  LxkLxu  "  LykLyu 


cpLtu  =  -g  +  LxkLxu  +  LykLyu  +  LzkLzu 
The  first  equation  becomes 


kLxu  =  u(0,y,z,t)  +  Lxg  +  LxcpLtu 


-  L"\..kL..u  -  L.-1L_kL_u 
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jvl 


Lxu  =  k_1[u(0,y,z,t)  +  L^g] 

+  k^L^cpI^u  -  k^L^Lyklyi 

-  k“\‘\  kL  u 

A  tm  4m 

u  =  u(0,y,z,t)  +  L"1{k"1[u(0,y,z,t)  +  L^g]} 

+  L"1k"1LycpUu  -  CVVY  kL  u 

A  A  V»  A  A  Jr  jr 

■  Lx1k'lLxlLzkLzu 

We  proceed  similarly  with  the  other  three  equations.  Then 

kLyu  ■  u(x,0,z,t)  +  L^g  +  L^cpL^  -  Ly^kl^u  -  Ly^kl^u 
LyU  =  k-1[u(x,0,z,t)  +  L^g]  +  k'^^cpL^  -  k'1Ly1LxkLxu 

-  k'VLzkLzu 

u  =  u(x,0,z,t)  +  L"1k_1[u(x,0,z,t)  +  L^g] 

+  Ly ^ k” 1 1”1  cpL^u  - 
The  third  equation  becomes 

klzu  =  u(x,y,0,t)  +  L^g  +  L^cpLj.u 
-  Ll1LvkLwu  -  L~^L  kl  u 

z  a  a  z  jr  y 

Lzu  *  k”1u(x,y,0,t)  +  k'^L’^g  +  k^L^cpL^u 

-  k’VLxV  -  k"VW 
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u  =  u(x,y,0,t)  +  L"1k"1u(x,y,0,t)  +  L“1k“1L“1g 

+  -  l;V'l-Vlxu 

-  L'VVV  kL,u 

€,  £.  J  J 


The  fourth  equation  is 


Ltu  =  -g/cp  +  (l/cp)(LxkLxu)  +  (l/cp)(LykLyu) 


+  l/cp(LzkLzu) 

u  =  u(x,y,z,0)  -  (g/cp)  +  L^(l/cp)LxkLxu 
+  L^d/cpJd  kL  u) 


Adding  and  dividing  by  four  we  get  u: 

u(0,y,z,t)  +  L~^  (k"^[u(0,y,z,t)  +  L^g]) 

+  L‘V1LvcpLtu  -  L"1k~1L"1L  kL  u 
x  xu  x  x  jr  y 

-  Lx'k''LxlLZkLzU 

+  u(x,0,z,t)  +  L"1k_1[u(x,0,z,t)  +  L”'g] 

+ 

+  u(x,y,0,t)  +  L”1k~1u(x,y,0,t)  +  L^VL^g 
+  ^lk'1L;1cpLtu  - 

-  LlVVV.kL.u 
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7)  Three-dimensional  linear  oaee:  Of  course  we  are  by  no  means  limited 
to  linear  equations;  the  decomposition  method  solves  highly  nonlinear 
equations.  However  we  consider  a  case  familiar  to  engineers  for  clarity 
and  simplicity.  Now  assuming  k,  c,  p  are  constants,  we  have  the  form: 

2  2 

a4v  u  -  ut  =  f 

Now  define  Lxx  =  32/9x2,  Lyy  =  92/9y2,  Lzz  =  92/9z2,  Lt  =  3/9t.  Then 
a  ^Lxx  +  Lyy  +  Lzz^u  Ltu  =  f 

Solving  for  each  linear  operator  term  in  turn,  we  have  the  four  equations 


a2L„u  = 

f 

+ 

L+u 

-  a2L  u 

-  a2L  u 

xx 

t 

yy 

zz 

a2L  u  = 

f 

+ 

L.u 

2 

-  a^L  u 

-  a2L  u 

yy 

t 

xx 

zz 

a2L  u  = 
zz 

f 

+ 

Ltu 

-  a2L  u 
xx 

-  a2L  u 
yy 

Lt“  ■  -f  *  +  aV  +  a2Lzzu 

or 

LXXU  ■  a'Zf  +  a'2|-tu  -  Lyyu  -  Lzzu 

V  ’  a'2f  +  “'V  -  LXXU  '  LzzU 

Lzzu  =  a‘2f  +  “'V  -  Lxxu  ■  V 

Lt“  -  -f  +  AJXU  ♦  Ayyli  ♦  a2Lzzu 


For  simplicity  let  us  allow  a  =  1  here.  Then 


L  u  =  f  +  L.u  -  L  u  -  L  u 
xx  t  yy  zz 


L  u  =  f  +  L..U  -  L  u  -  L  u 
yy  t  xx  zz 


mm 


ml 


L  u  =  f  +  L.u  -  L  u  -  L  u 
22  t  xx  yy 


L.u  =  -f  +  L  u-L  u-L  u 
t  xx  yy  22 


Operating  with  the  inverse  operators  (two-fold  integration  for  x,  y,  2 
and  single  for  t), 

U(x,y,2,t)  "  A  ♦  Bx  +  L;Jf  +  CL^Lt  -  L^Lyy  -  L^LJu 
u(x,y,z,t)  ■  C  +  Dy  ♦  L“’f  +  [LyyLt  -  L'Jlxx  -  L;Jl„]u 
u(x.y,2,t)  -  E  ♦  Fz  +  L'Jf  +  [i;’Lt  -  L']lxx  -  L-JLyy]u 
u(x.y.z.t)  =  G  -  l-V  +  [L;1^  ♦  L'\y  *  I^L^U 


Again,  adding  these  equations  and  dividing  by  four,  we  have  u  on  the 

00 

left  side  and  a  complex  expression  on  the  right  side.  Letting  u  =  Yu 

n=0  n 

we  identify  uQ  as  including  terms  involving  the  boundary  conditions 
and  the  forcing  function  f.  A,  B,  ...,  G  are  chosen  to  satisfy  boundary 
conditions.  Now  we  have 


uQ  =  (1/4){A  +Bx+C+Dy+E+F2+G 

-  +  lyy  +  L«  -  L^f} 

un+l  ”  ^xx^t  "  ^xxLyy  “  *'xxL22  +  Lyy  **t  “  ^yyLxx  ”  *"yyLzz 

+  L22Lt  "  L22Lxx  "  L22Lyy  +  Lt  Lxx  +  Lt  Lyy  +  Lt  L2z^un 

n-1 

for  n  >  0.  The  expression  4>n  =  £  u..  is  the  approximation  to  u,  i 
the  correct  solution  is  the  sum  of  all  the  components  u..  but  since  the 
series  converges  so  rapidly  that  a  few  terms  are  sufficient  for  a  good 
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approximation,  some  n  terms  will  be  enough  or  <|>n  =  Ug  +  u-j  +  •••  +  un~^ 
n-1 

or  l  u.  will  serve  as  an  approximation  to  u.  Easily  computable  and 
i=0  1 

accurate  solutions  have  been  obtained  with  small  n  for  differential  and 
partial  differential  equations,  even  when  nonlinear  terms  or  stochastic 
processes  are  included.  Numerical  computations  clearly  demonstrate  the 
convergence  to  desired  accuracy. 


8)  Tuo-dimeneional  homogeneous  heat  equation:  As  an  illustration  of 
procedure,  we  consider  the  specific  example 

V  U  -  Ut  =  0 

p 

with  given  conditions  u(0,t)  =  t,  u(x,0)  =  x  /2,  and  3u/3x[x=Q  =  0; 
we  now  have 

u  =  u(0,t)  +  LxjLtu 

u  =  u(x,0)  +  L^Lxxu 

Adding  and  dividing  by  two,  we  have 

oo 

u  =  (l/2)[u(0,t)  +  u  (x  ,0)  ]  +  (l/2)[L^Lt  +  L'1^]  £  un 


where  we  have  written  u  *  Y  u  . 

n 

n=u 


uQ  =  (l/2)[t  +  (x2/2)] 
ui  ■  <,/2)[L«Lt  *  l;’l 


XX 


Now  we  can  identify 

]u0 


>.<  f.t  (.-•  *.'  ( , 


,»•  .4,  . i.  i 


*'»•!  VI  '.'I  *♦!  Vl’t'l'i*!  '•*»  t*fl  »*#  i*»  t'l  ill 


u2  -  (i«ici^it  * 


un+l  *  ^2^LxxLt  +  Lt  Lxx^un 


n  >  0 


This  specific  case  with  the  given  conditions  has  been  chosen  so  that 


[L'jL,  +  Li1 L  ]  is  an  identity  operator,  a  special  case  making  the 

XX  t  1  XX 


calculation  particularly  simple  for  the  sake  of  clarity.  Thus  we  see 


[LxxLt  +  LtlLxx][u(0»t)  +  u(x*0^ 


-  L;iLtu<°*t>  *  +  LilLxx“(°'t)  +  Li,Lxxu'x-°> 


-  L«Lt<t>  +  L«Lt<x2/2>  +  ^\x<‘>  ♦  L-1lxx(x2/2) 


=  (x2/2)  +  t  =  u(x,0)  +  u(0,t) 


In  the  general  case  one  has  to  compute  the  effects  of  these  operators 


repeatedly,  but  it  is  still  straightforward.  Now  for  the  case  being 


considered 


uQ  =  (1/2)  [t  +  (xV2)] 


U1  ■  <^2)[l«Lt  *  LilLxx^u0  ‘  '1/4)[t  +  (x2/2,] 


u2  -  (l/2)[LxJl.t  +  L;,Lxx]u,  =  (l/8)[t  +  (x2/2)] 


u„  =  (l/2ntl)[t  +  (x2/2)] 


u  -  [t  *  (x2/2)]  l  (l/2n+1 ) 
n=0 


is  the  solution  obtained  from  summing  the  un  from  0  to  °°.  Denoting 


m 


mfo i 


i 


i!l$§ 


m 


n 

mm 

m 

a 


Us&SA 


m 


$$g 

K'*  Hi 


m 

m 


urn 


livivi  i' 


mm 

OT.I 


m 


id 


mm 


m 


54 


jctjwjw  MjrsunCTTonwrwvwvro* 


an  n-term  approximation  by  <j>n, 

4>1  =  0.5(t  +  x2/2)  =  uQ 

<f>2  =  0. 75{t  +  x2/2)  =  uq  +  u-j 

(j>3  =  0 . 87  5  ( t  +  x2/ 2)  =  Uq  +  u-j  +  u2 

u  *  Tim  <j>n  =  u-|  +  u-|  +•••+  un-1  =  (l)(t  +  x2/2) 


2 

which  shows  the  improving  approximations  to  the  correct  solution  t  +  x  /2 
as  n  increases.  We  see  the  convergence  of  the  coefficients  .5,  0.75, 
0.875...  to  1.0  illustrating  the  improving  approximation.  Solutions 
using  the  decomposition  method  are  analytic  but  not  closed  form  -  the 
closed  form  in  this  case  arose  from  a  problem  chosen  for  the  purpose  of 
clear  illustration.  Note  that  six  terms  yield  the  solution  to  better  than 
98%  and  with  10  terms,  the  approximation  is  within  99.9%  of  the  correct 
value. 

Addition  of  stochasticity  in  forcing  functions  or  coefficients  results 
in  the  series  <J>n  becoming  a  stochastic  series  from  which  statistics  are 
easily  obtained  by  averaging. 


Example:  Heat  conduction  in  a  uniform  beam:  For  one-dimensional  flow 
parallel  to  the  x-axis,  the  equation  is  u.  =  ku  where  k  is  the 

t  X/\ 

thermal  diffusivity.  Taking  k  =  1,  we  have  u.  =  u  if  heat  is 

v  X  A 

transferred  only  by  conduction,  and  no  internal  source  is  present. 

Boundary  conditions  describing  the  thermal  conditions  on  the  surface 
of  the  solid  and  the  initial  temperature  distribution  are  necessary  along 


ft# 


with  the  heat  equation  to  determine  the  temperature  distribution  u(x,t). 
(For  a  three-dimensional  solid  with  an  insulated  end  at  x  =  0,  the 
condition  on  u  is  u  (0,y,z,t)  =  0  because  the  heat  flow  across  the 

A 

surface  would  be  proportional  to  ux(0,y,z,t).) 

Suppose  we  have  a  bar  whose  temperature  distribution  is  described 
by  u(x,t).  Let  the  bar  be  initially  at  temperature  f(x)  =  sin  x  and 
assume  the  ends  x  =  0  and  x  =  tt  are  kept  at  temperature  zero.  We 
have  therefore  the  conditions  u(x,0)  =  sin  x,  and  u(0,t)  =  uU.t)  =  0. 
Writing  Lxxu  *  L^u  and  applying  the  inverse  operator  L~x  to  both  sides, 
we  get 


u  =  A  +  Bx  +  L-xLtu. 


m. 

i 


$$ 


p 

m 


A  and  B  must  be  zero  to  satisfy  the  given  boundary  conditions  at  the 
ends  so  that  uQ  is  zero  in  this  equation.  Whenever  this  happens,  that 
particular  equation  does  not  contribute  to  the  solution.  (The  decomposition 


method 


requires  a  non-zero  uQ.)  Now  applying  L“‘  to  Ltu  =  Lxxu 


and  satisfying  the  initial  condition,  we  have 


u  =  u(x,0)  +  Lj;  Lxxu 


so  that  Ug  =  u(x,0)  =  sin  x.  We  have 


“  *  u0  +  L*x  J0  un 


U1  ’  Lt'Lxxu0  ■  Lt’Lxx  s("  x  =  sin  x 

u2  *  LtlLxxul  ’  Lt1'-xx(-t  5in  x)  '  (t2/2)  si"  x 
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2  -t 

The  solution  is  u  =  (1  -  t  +  t  /2  -  ...)  sin  x  =  e  sin  x  which  is 
easily  verified. 

It  is  interesting,  as  a  simple  test,  also  to  consider  the  trivial 
example  of  steady-state  temperature  in  a  bar  with  ends  at  x  =  0  and 
x  =  1,  maintained  at  u  =  0  and  u  =  T  respectively.  Then  u(0)  =  0 
and  u ( 1 )  =  T.  Since  we  are  considering  steady  state,  we  have  no  time 
derivative  and  only  uxx  =  0  so  that  u  =  Ax  +  B  is  the  complete 
solution  instead  of  just  the  first  component  Uq.  Applying  the  conditions, 
i.e.,  the  known  temperatures  at  the  ends,  we  have  B  =  0  and  A  =  T  so 
that  u  =  Tx  is  the  well-known  solution. 


9)  Conductivity  dependent  on  x:  Consider  a  rod  of  length  l  lying 
along  the  x  axis  with  its  left  end  at  the  origin.  Suppose  its  density 
is  a  function  of  x,  its  cross-section  is  A,  its  mass  per  unit  length 
is  p,  specific  heat  is  c,  and  conductivity  is  k-j .  Let  u(x,t) 
represent  temperature  at  a  point  x  at  time  t  and  let  F(x,t)  be  the 
amount  of  heat  per  unit  cross-section  per  unit  time  passing  x  to  the 
right.  Then,  the  cylinder  of  base  area  A  between  the  points  x  and 
x  +  6  has  heat  supplied  to  it  at  the  rate  cpSAu^  =  AF(x,t)  -  AF(x  +  6,t)  = 
-AFx6  where  higher  powers  of  6  are  dropped  since  6  is  small.  Since 
F  =  -k,u  in  heat  conduction,  F  =  -k,u  .  Now  we  have  the  model  equation 

I  «  ^  i  Aa 

<c»/kl>ut  =  \x 

or  letting  k  =  cp/k-j,  we  have  the  heat  equation 
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where  conductivity  k  is  not  a  constant  but  is  dependent  on  position. 
Suppose  we  have  initial  boundary  conditions  u(x,0)  =  g(x)  for  0  <  x  <  1 
and  u(0,t)  =  h-j(t),  u(l,t)  =  h2(t)  for  t  >  0.  Then 

uQ  =  (l/2)[g(t)  +  h-j (t)  +  x(h2(t)  -  H1  (t) )3 

from  which  the  following  terms  can  be  computed  analogous  to  the  previous 
procedure  when  k  was  assumed  constant.  However  we  now  have 


Lxx“  ■  k<x)Ltu 


u  =  A  +  Bx  +  L^k(x)Ltu, 


k(x)Ltu  =  Lxxu 


Ltu  =  k~  (x)Lxxu 
u  =  C  +  L^k”1  (x)Lxxu, 


Uq  =  A  +  Bx 


uQ  =  c 


so  that  if  neither  uQ  vanishes 

u  =  (1/2)[A  +  Bx  +  C]  +  (l/2)[L"j|k(x)Lt  +  L'V1  (x)Lxx]u 
uQ  =  (1/2)[A  +  Bx  +  C] 

V,  ■  0/2)[L^(x)Lt  ♦  L-tH-'(x)Lxx]an 


Now  we  can  consider  various  specific  examples.  Our  only  objective  here  is 
to  demonstrate  solutions  of  problems  also  solvable  by  other  methods  for 
easy  verification  and  to  show  the  method  clearly  and  the  simplicity  of 
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the  computation  involved.  Let  h1  =  h2  =  0  and  g(x)  =  sin  ttx.  Then 
one  equation  is  sufficient.  If  k  =  1, 

uQ  =  sin  ttx 

-1  2 
u-j  =  Lt  Lxxsin  ttx  =  -ir  t  sin  ttx 

u2  =  *-xx ( -TT^t  sin  ttx)  =  (Tr^t^/2)  sin  ttx 

u  =  [1  -  (ir2t)  +  (ttV/2)  -  ••*]  sin  ttx 

At  this  point  we  happen  to  be  able  to  identify  the  bracketed  series  as  e 
so  we  can  write 

— tt  2 1 

u  =  e  sin  ttx 

Identification  of  a  closed  form  is  not  necessary.  We  have  an  analytic 
approximation  that  can  be  calculated  numerically  to  see  it  converges  to 
the  correct  result.  Thus  the  temperature  u(x,t)  can  be  calculated 
exactly  at  any  x  and  t.  To  consider  more  general  cases  we  can  allow 
a  given  k(x),  e.g.,  k(x)  =  1  +  x,  and  compute  a  series  for  u(x,t). 
The  method  can  now  be  generalized  to  parabolic  equations  such  as 

V  u  -  kUj.  =  g 

where  k  =  k(x,y)  or  k(x,y,z)  and  g  =  g(x,y,z,t)  or  even  stochastic 


cases. 
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10)  nonlinear  heat  equation  in  one-dimension: 


(3/3x)[k(u)3u/3x]  =  cp(3u/3t) 


u(0,t)  =  u 


1 


u(x,0)  =  u2 


Here,  we  can  let  k(u)  be  any  function  -  a  linear  function  kg  +  k-jU,  or 


a  +  Bu  +  yu  »  or  sinh  u.  Suppose  we  consider  the  first.  We  have 


Lxk(u)Lxu  =  cpltu 


or 


LxkoLxu  +  LxkluLxu  •  «Ltu 


The  procedure  is  straightforward.  We  have  two  linear  terms  ko^x^x,J 
and  cpLtu  and  the  nonlinear  term  k-|Lx(uLxu).  In  the  nonlinear  term, 


we 


replace  uL  u*=uu  by  l  A  where  the  A  are  generated  for 

A  A  _ _ n  II  "I 

n=u 


Nu  =  uux<  We  then  solve  for  the  two  linear  terms  in  turn  using  the 


inverse  operators  L~^  and  .  We  define  Uq  as  usual  and  assume 


the  decomposition  u  =  £  u  .  The  components  u, ,  u9, 

n=0  n  1  L 


now  are 


determinable. 

Solving  for  the  linear  term  L  knL  u  and  substituting  a  for 

X  v/  X 


cp, 


k0Lxu  ■  u(0't)  +  Lx1<lLtu  •  L^LxkluLxu 


Lxu  -  k-1u(0.t)  +  koVaLtu  -  koVLxkluLxu 


M 


WMM? 


m 


u  =  u(O.t)  *  L-Jk^(O.t)  ♦  L->o'L^aLtU 
Lxxk0  LxxLxxkluLxxu 

u  -  «(0,t)  ♦  L^Jk-’utO.t)  *  k0’al^L^Ltu 
kxxk0  kxxkxxklukxxu 


Now,  solving  for  the  remaining  linear  term 


Ltu  =  a'lLxxk0Lxxu  +  “'lLxxk1Lxxu 
u  =  u(x,0)  *  o'1k0L;,LxxL)(XU  ♦  a'Vt’LxxuLxxu 


Adding  the  two  equations  for  u  and  dividing  by  two: 

Uq  =  (l/2){u(0,t)  +  u(x,0)  +  L’Jk^utO.t)} 

u "1  “  (1/2^{k0laLxxLxxLtu0  "  Lxxk0lLxxLxxklu0Lxxu0 

+  a  1 k0LtlLxxLxxu0  +  a  lklLtlLxxu0LxxU0} 

u2  =  (1/2^{k0laLxxLxxLtul  '  LxxkOlLxxLxxk(ulLxxuO  +  u0Lxxul) 

+  «‘V;\xLxxul  +  a"1klLilLxx(u0Lxxul  +  ulLxxu0)} 
so  that  components  are  determinable  as  necessary. 
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11)  Non-uniform  Conduction  or  Nonlinear  Boundary  Conditions:  Suppose  the 
heat  conducting  material  is  not  uniform  so  that  the  coefficients  in  the 
equation  are  discontinuous  at  one  or  even  several  points.  We  can  then 


divide  the  total  interval  [0,2,]  at  the  points  of  discontinuity 
into  subintervals.  If  the  temperature  and  the  heat  flow  are  continuous 
at  the  points  e^,  we  have 

u(ei  -  0,t)  =  u(ei  +  0,t) 

k(e1-  -  0)3u/8x(Ei  -  0,t)  =  k(e.  -  0)3u/3x(e.  +  0,t) 

If  heat  is  being  radiated  at  the  x  =  0  cross  section  of  the 
conducting  material  with  the  temperature  T(t),  we  have  the  nonlinear 
boundary  condition: 

k  3u/3x(0,t)  =  o[u4(0,t)  -  T4(0,t)] 

Such  nonlinear  boundary  conditions  can  also  be  handled  by  the  decomposition 
method.  The  procedure  is  simple.  We  use  the  nonlinear  boundary  equation 
given  just  like  a  differential  equation  again  using  the  approximation 
<J;n  -  for  example,  the  three-term  approximation  4>3  -  as  the  solution  of 
the  boundary  equation  which  results  in  evaluation  of  the  constants  of 
integration  in  <j>3.  A  large  number  of  examples  of  this  type  have  now 
been  verified  for  accuracy. 

The  power  and  convenience  of  the  decomposition  method  is  shown  by  the 
following  examples  for  comparisons  with  usual  procedures. 

2 

Consider  the  homogeneous  heat  conduction  equation  ut  =  a  uxx  with 
u(x,0)  =  <j>(x)  and  u(0,t)  *  u(2,t)  =  0. 


ill 


Applying  the  well-known  'separation  of  variables  technique 
assuming  <j>(x)  is  continuous,  bounded,  possesses  piecewise  continuous 
derivatives  and  satisfies  <j>(0)  =  <$>(£)  B  0.  (Then  u(x,t)  will  be 
continuous  for  t  >0.)  Let  u(x,t)  =  X(x)T(t).  Then  we  have 

(1/a2)  (T*/T)  =  (XVX)  =  -X  =  constant 


which  leads  to  the  equations 


X"  +  XX  -  0 


An  •  (irn/S.)‘ 


Xn(x)  =  sin(irnx/£)Z 
Tn(t)  -  Cne-aV 


where  the  Cn  are  constants.  Thus 


X(0)  =  X(£)  =  0 


where  n  =  1 ,2,3,, 


u(x,t)  =  l  X  (x)T  (t)  =  l  c  e"a  Xntsin(un/£)x 
n=0  n  n  n=0  n 

satisfies  the  homogeneous  boundary  conditions.  To  satisfy  the  initial 
conditions 


<|>(x)  =  u(x,0)  -  l  C  sin(irnx/£) 
n=l  n 

which  is  the  Fourier  sine  series  so  that  the  c  are  Fourier  coefficients 

n 

which  we  can  evaluate 


=  =  2/S,  <J>(e)si  n(Trne/S,)de 

n  n  jn 


[It  can  be  shown  that  the  series  converges  and  is  appropriately  (term-wise 
twice)  differentiable  with  respect  to  x,  and  satisfies  the  differential 
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equation  in  0  <  x  <  £,  t  >  0.]  We  will  see  below  that  solution  by 
decomposition  is  very  much  simpler  and  that  it  results  in  a  much  more 
easily  computed  result.  Before  we  proceed  with  decomposition  we  would 
like  to  comment  that  this  solution  can  also  be  written  as 


u(x,t)  =  G(x,e,t)<j>{c)d£ 

Jn 


when  the  Green's  function  G  is  given  by 


G  =  (2/£)  l  exp{-(7rn/£raZt}  sin(7mx/£)  sin(7rne/£) 
n=l 

which  can  be  seen  by  looking  at  u(x,t),  substituting  values  of  cn- 

Since  the  series  for  converges  uniformly  for  t  >  0  with  respect 

to  e,  the  summation  and  the  integration  can  be  interchanged. 

2 

If  instead  of  ut  =  a  uxx>  we  consider 

ut  =  a2uxx  +  f(x,t) 

i.e.,  a  more  general  equation  with  a  forcing  function  and  assume  given 
conditions  as: 


u(x,0)  =  u(0,t)  =  u(£,t)  *  0 


the  solution  becomes 

(t  fl 

u(x,t)  =  G(x,£,t-x)f(£,x)d£dx 

Jo  Jo 

G  =  (2/£)  l  e“(7Tn^)  a  (t-T)sin(Trnx/£)sin(TTne/£) . 
n=l 

We  can,  if  we  wish,  apply  decomposition  to  the  separated  equation, 
but  we  emphasize  that  there  is  no  need  for  the  separation  anymore  since 
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we  can  do  the  partial  differential  equation  directly  as  demonstrated  by 
our  examples  so  far.  The  decomposition  result  is  very  much  simpler  and 
much  easier  to  compute  than  the  result  above.  If  we  add  nonlinear  or 
stochastic  terms  as  well,  then  the  decomposition  method  works  as  easily 
and  yields  a  much  superior  result  because  no  artificial  simplifications 
are  forced  on  us. 


12)  Self- shadowing  Latticework  Structures  in  Orbit:  An  orbiting  latticework 
structure  -  a  truss  -  undergoes  time-dependent  heating  by  the  sun  and  the 
earth,  as  well  as  by  power  systems,  heat-exchangers,  and  electronics  on 
the  space  station.  Meanwhile  heat  is  being  radiated  by  the  truss  itself. 

The  temperature  distribution  T(a,t)  along  a  truss  member  is  obtained 
by  solving  conservation-of-energy  equations  -  ordinarily  accomplished  by 
finite-element  methods. 

Since  the  specific  equation  is  not  significant  in  this  discussion  and 
the  heating  history  depends  on  the  orbit,  the  shadowing  by  the  rest  of  the 
structure,  etc.,  let  us  first  consider  a  simple  case  with  an  isothermal 
truss  member.  The  truss-member  has  temperature  T(t)  which  must  be 
found  from 

3T/3t  +  kT4  =  g 


where  k  depends  on  density,  emissivity,  volume,  emitted  radiation 
surface  area,  and  the  Stefan-Boltzmann  constant  and  g  depends  on 
Incident  radiation  surface  area  and  heating  rate.  The  equation  can  be 


rewritten  In  our  form 


LfT  =  g  -  kr 


Applying  the  inverse  operator  as  usual,  we  have 


T  =  T( 0)  +  L^g  -  kl^M 


and  by  decomposition 


T„  -  T(0)  +  L-’g 


T,  *  -kL-’A^T4) 

T2  =  -kL^A^T4} 


Vi  ’  -kLi’A„(T4} 


where  l  T  ■  T.  An  n-term  approximation  will  then  be  given  by  4>  = 
n=0  i 

or  Tq  +  T'i  +  •  •  •  +  Tn  -j  which  we  have  shown  how  to  calculate  above  once 

4 

we  know  the  An>  The  An  polynomials  for  the  quadratic  nonlinearity  T 
are  available  from  the  published  literature*  and  are  given  by: 

A  =  T4 

rt0  'o 


A1  =  4T0T1 

A2  =  4T^T2  +  6TqT^ 

A3  =  4TqT3  +  4T^Tq  +  1 21^1^0 


See  reference  [1,2,23]. 


tew 

Um 

mm 

m 


66 


r 

fUVvVaVVVVV 


A4  =  +  4TqT4  +  6T2T2  +  12TqT.|T3  +  12T2TQT2 


A6  .  4TpTg  +  4T?T2  ♦  12t|TiT4  ♦  lZT^Tj  ♦  IZT^Tj  f  lZT^T, 


Thus  the  system  is  completely  calculable. 


When  we  have  a  more  complex  situation  with  multiple  truss-members, 


systems  of  coupled  equations  arise.  These  are  solved  with  equal  facility 


by  decomposition. 


Also,  since  heating  will  not  be  uniform,  the  isothermal  approximation 


is  obviously  poor  and  we  will  have  nonlinear  partial  differential  equations 


such  as 


9T/3t  -  «92T/9x2  +  ST4  =  g 


where  a  depends  on  thermal  conductivity,  cross-sectional  area,  density, 
specific  heat,  and  volume  and  B  depends  on  emissivity,  emitted  radiation 


surface  area,  and  the  Stefan-Boltzman  constant.  In  order  to  solve  for 


T  we  write 


LtT  -  “LxxT  *  ST  =  g 


This  we  now  solve  by  decomposition  by  first  solving  for  the  two  linear 


operator  terms;  thus  we  get  the  two  equations 


LtT  "  9  +  “LxxT  -  «T 


“LxxT  =  -9  +  LtT  ♦  ST 


From  equation  (1)  if  a,  B  are  constant 


T  =  T(0)  ♦  L-’g  *  aL-’L^T  -  BL^T4 


"Self-Shadowing  Effects  on  Thermal  Structural  Response  of  Orbiting 
Stre-ses,"  «7.  Mahoney  and  E.  Thornton ,  J.  Spacecraft,  24,  no.  4,  pp.  342-348. 
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From  equation  (2) 


T  e  A  +  Bx  -  a-’L-Jg  +  c-'l'JlJ  ♦  cf'el^T4 


(4) 


With  the  decomposition  of  T  into  components  to  be  found,  i.e., 

4 


T  =  l  Tn  and  representing  the  nonlinear  term  V  using  the  A 
n=0  _  n 


4  A  _  A 

polynomials  for  T  ,  i.e.,  T  =  £  A  {T  }  we  have 


n=0 


A  = 

rt0  'o 


A,  =  4TJT, 


A2  =  4T2T2  +  6T2T2 


A3  '  4T0T3  +  4T?T0  +  ,2TiTlT2 


A4  =  T?  +  4T0T4  +  6T0T1  +  ,2T0T1T3  +  121T0T2 


A5  =  4T0T5  +  4T?T2  +  ,2T0T1T4  +  12T0T2T3  +  12T1T0T3  *  ,2T2T0T 


From  equation  (3) 

T0  -  T(0)  +  L-’g 

Tn+1  =  aLt  LxxTn 
From  equation  (4) 


A 


n 


Tq  =  A  +  Bx 

Tn+1  ”  "a  Lxx9  +  a  LxxLtTn  +  a  B^xxAn 
Here  A  and  B  are  determined  by  allowing  Tq  or  ,  the  first  or 
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one-term  approximation  to  T,  to  satisfy  given  conditions  T(0,t)  and 

T(t,t)  for  the  truss  member  temperature.  Since  TQ  is  now  known  in 

both  equations,  the  one-term  approximate  for  T  is  ,  which  is  half 

* 

the  sum  of  the  two  4>-j '  s 

Now  we  can  calculate  the  two  u^'s  and  continue  in  the  same  manner 
until  we  have  an  n-term  approximation  <f>n,  i.e.,  an  approximation  to 

the  correct  temperature  T  which  becomes  more  accurate  as  we  go  to  more 
terms  until  there  is  no  further  change. 

It  will  now  be  clear  that  the  decomposition  method  is  not  only 
applicable  to  generic  space  structure  thermal  problems  but  provides  an 
extremely  useful  straightforward  method  as  compared  with  ordinary  numerical 
procedures. 

Consider  some  examples  carried  out  to  a  complete  solution  to 
demonstrate  the  validity  of  these  procedures. 

EXAMPLE  1: 

Consider  the  application  of  decomposition  to  an  equation 

u  -  u  =  0  modelling  a  square  plate  with  0  <_x  <_  tt/2  and  0  <  y  <  n/2 
xx  yy 

with  the  given  conditions: 

u(0,y)  «  0  u(ir/2,y)  =  sin  y 

u(x,0)  =  0  u(x,u/2)  =  sin  x 


In  some  special  cases,  one  equation  may  not  contribute  because  the  1 
term  is  zero.  Then  we  simply  work  with  the  remaining  equation  and  the  8 
adding  and  dividing  is  unnecessary.  See  example  (4)  on  a  following  page.  ij 


4  »'  -.1  1 


Let  L  -  3  /3x  and  L  =  3  /3y  and  write  the  above  equation 
xx  yy 

as  L  u  =  L  vu.  As  usual  in  the  decomposition  method,  we  solve  for 
xx  yy 

each  linear  operator  term,  L  u  and  L  u,  in  turn  and  then  apply 

xx  yy 

the  appropriate  inverse  to  each. 

L;xLxxu  *  u  -  clklW  '  c2k2(),)x  '  LxxV 


'xx  yy 


LyyV  '  u  '  c3k3<x>  -  c4k4(x>*  -  LWLXXU 


u  =  clkl(>>  +  c2k2<*>x  +  LxxLyyu 

u  =  c3k3(x)  +  c4k4(x)y  +  LyyLxxu 

Define  <}>x  =  c^^y)  +  c2k2(y)x  and  <|>  =  c3k3(x)  +  c4k4(x)y 
to  rewrite  (1)  and  (2)  as 


u  =  <b  +  L~  L  u  (3) 

yx  xx  yy  ' 

u  =  4>  +  L”^L  U  (4) 

y  yy  xx  ' 

One-term  approximants  to  the  solution  u  are  uQ  =  <f>x  in  (3)  and 

u0  =  ^y  Tw0-term  approximants  are  Ug  +  u^  where  u-j  =  L~xLyyu 

in  (3)  and  L^LxxUg  in  (4),  etc.  Thus  un+1  =  LxxLyyun  in  (3)  and 
LyyLxxun  1n  (4>  for  n  i  °- 

For  the  x  conditions  u(x, 0)  =  0  and  u(x,tt/2)  =  sin  x  applied 
to  the  one-term  approximant  uQ  =  c-j k-j (y )  +  c2k2(y)x,  we  have 

Cjk^y)  =  0 

c2i ^(yW2  ~ sin  y 


m 

mu 


n**ViV 

IVMi't 
W'i  't 

Eh&M'i 

5# 

BSP 


m 

ii 


m 

P 

Wy'i 


■MS' 


SfflS 

PVV* 


I 


or  c2  =  2/tt  and  k2(y)  =  sin  y. 

For  the  y  conditions  u(x,0)  =  0  and  u(x,tt/2)  =  sin  x  applied  to 
u0  *  c3k3(x)  +  c4k4(x)y,  we  get 

c3k3(x)  =  0 
c4k4(x)Ti/2  =  sin  x. 

Thus  c4  =  2/ir  and  k4(x)  =  sin  x. 

★ 

If  a  one-term  approximant  were  sufficient,  the  solution  would  be 
4>i  =  (1/2) { (2/tt)x  sin  y  +  (2/ir)y  sin  x) 

The  next  terms  for  (3)  and  (4)  respectively  are 

U1  ■  L«Vo  =  L«LyyCc2x  s1n  ^ 

U1  ■  LyyLXxu0  ’  LyyLx*[V  si"  XJ 
We  continue  to  obtain  u2,u3>...  .  Clearly,  for  any  n, 

un  =  (L™Lyy>°u0  ’  c2(si"  y > (-1  >nx2n+1 /{2n+1 } : 

u„  -  (l;Jlxx)"u0  =  c4(sin  x)(-l)nyZn+1/{2n+1>: 

Letting  <j>m  represent  the  m-term  approximant,  we  have  for  the  two  cases: 


sin  y  V  (-l)nx2n+1/(2n+l): 

(5) 

n=0 

sin  x  (-l)ny2f1+1/(2n+l): 

(6) 

We  can  now  apply  the  conditions  4 >  (n/2,y)  =  sin  y  for  (5);  thus 


An  n-term  approximant  is  J  u. 

n  i=0 
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c1k1 (y)  -  0 


c2  sin  y  ^  (v/2)Zn+'/(2n+l )!  =  s 


EXAMPLE  II: 


If  we  modify  our  example  to  u  -  u  +  f(u)  =  0  we  let 

xx  yy 

oo 

f(u)  =  1  A  (f(u)},  i.e.,  f(u)  is  represented  by  the  polynomials 

n=0  n  n 

generated  specifically  for  f(u).  Then  the  equations  (3)  and  (4)  in  the 
previous  example  become 


and 


U  ’  *x  +  LxxLyy  j0  un  '  jQ  Vf<u» 


u  =  <f>  +  L_1L  y  u  +  L"1  y  A  (f (u) } 
y  yy  xx  ^  n  yy  n£0  V 

and  computation  proceeds  as  before. 


-1 


EXAMPLE  III: 


Now  consider  u 

A  A 

u(x,0)  =  sin  x 
ut(x,0)  =  0 


utt  on  0  £  x  £  it  with  t  £  0  specifying 
u(0,t)  =  0 
u(it,t)  =  0 


Decomposition  yields 

u  «  Ci k-j (t)  +  c2k2(t)x  +  L^Ltt  l  up  (1) 

n=u 

u  -  c3k3(x)  +  c4k4(x)t  ♦  l;{lxx  l 

n-u 

In  (1)  the  Uq  =  c^k^(t)  +  c2k2(t)x  term  vanishes  when  the  boundary 
conditions  u(0,t)  =  0  and  u ( tt , t )  =  0  are  applied.  Hence  in  this 


(2) 


case  (1)  does  not  contribute  and  we  need  consider  only  (2).  Here 

u0  =  c3k3^x^  +  c4k4^x^  and  applying  the  conditions  u(x,0)  =  sin  x  and 

ut(x,0)  =  0,  we  have 


Uq  =  sin  x 


ui  =  LnLnxuo  =  <-t2/2>  s1"  x 

u2  *  LttLxxui  ‘  (t4/4l)  sin  x 


u  =  sin  x  cos  t 


Here,  we  recognized  the  series  for  cos  t  but  ordinarily  use  the 
series  and  see,  most  effectively,  in  numerical  solutions  that  the 

n-1 

solutions  stabilize  as  n  increases,  i.e.,  the  approximant  =  l  u. 

co  n  i=0  1 

converges  to  the  correct  solution  u  =  £  u„. 

n=0  " 


EXAMPLE  TV: 


We  point  out  an  example  where  one  linear  operator  term  does  not 
contribute.  Consider  the  example  u  =  uf  with  given  conditions 

u(x,0)  =  sin  ttx  for  0  <  x  <  1 


u (0 , t )  =  u(l,t)  =0  for  t  >  0 


We  write 


Applying  the  L~]  to  (1)  we  have 


u  =  u0  +  LxxLtu 


Uq  =  A(t)  +  B(t)x  =  0 


because  of  the  conditions  u(0,t)  =  u(l,t)  =  0;  hence  this  equation 
doesn't  contribute  and  under  the  rules  of  the  decomposition  method  is  not 
considered  in  determining  u^. 

From  (2)  applying  L”1* ,  we  have 


u  =  u0  +  Lt  Lxxu 


un  =  u(x,0)  =  sin  ttx 


U1  ■  Li'Lxxu0  =  Ltllxx  si"  *x 
2 

=  -tt  t  sin  ttx 


u2  "  L't  Lxxul 


Summarizing  the  series,  we  determine 

— TT 2  t  . 

u  =  e  sin  ttx 


is  the  exact  solution  since 


u  *  Uq  +  u^  +  *•* 


=  sin  ttx  -  it  t  sin  ttx  + 


=  [1  -  tt  t  +  •••]  sin  ttx 


where  the  bracketed  quantity  is  the  series  for  e* 
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If  we  changed  the  x  interval  to  (0,tt/2)  and  specified 
u(x,0)  a  sin  x,  u(0,t)  =  0,  and  u(n/2,t)  =  e-t,  we  get  a  case  which 
is  less  convenient  computationally  because  of  the  integration  "constants." 

We  no  longer  have  un  =  0  for  the  L  u  equation  and  must  consider 
both  equations  in  determining  the  initial  term  Uq. 

Now  Uq  =  (2/ir)xe_t.  When  we  determine  u,,  we  get 
-(2/Tr)(x^/3!)e_t  +  c^(t)x  +  C2U)  because  of  the  two-fold  integration 
represented  by  L”x. 

When  x  =  0,  u  =  0;  hence  C2  =  0.  When  x  =  n/2,  u  =  e  S 
hence  (2/tt)(-jt/2)3(1/3!  )e_t  +  c^ (tt/2 )  =  e-t  determines  Cp  etc.,  and 
we  see  this  case  will  become  numerically  more  complex.  The  solution 
u  -  e-^  sin  x  is  still  valid  but  not  as  simple  to  obtain. 

We  might  also  note  that  if  u(x,0)  is  specified,  then  uxx  is 
specified  at  t  =  0  and  therefore  ut(x,0)  is  specified  through  the 
equation  so  we  cannot  independently  specify  ut(x,0).  Thus  the  conditions 
cannot  be  assigned  arbitrarily  and  must  of  course  be  physically  appropriate. 
We  can  specify  here  the  temperature  u  or  the  heat  flow  (derivative  ux) 
at  both  ends  or,  say,  temperature  at  x  =  0  and  heat  flow  at  x  =  1  (rod 
with  uninsulated  end). 


13)  Transient  Heat  Transfer:  K.  N.  Shukla  and  L.  Mani  (AIAAJ ,  October  1984) 
considered  thermal  constriction  resistance  with  arbitrary  heating  in  a 
convectively  cooled  plate.  Precise  heat  transfer  calculations  are,  of 
course,  important  to  achieve  reliability  in  space  structures.  The  two- 
dimensional  model  considered  in  the  above  reference  can  be  written  in 
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our  notation  as:  Lxxe  +  L  6  =  kL^e  where  txx  =  3^/3x^,  lyy  =  d^/dy^ . 

We  write  1  (mediately 

Lxx«  '  kLte  -  V 

V  *  kLt6  -  Lxxe 
Lte  -  k-\xe  ♦  k-’Lyye 

With  our  inversions 

6  .  A  +  Bx  +  kL^S  -  L"xLyye 

e  -  c  +  Oy  +  kL;jLte  -  L;’Lxxe 

e  »  e(t  -  o)  +  ♦  k-’i^e 

where  the  functions  A,  B,  C,  D  are  determined  from  specified  boundary 
conditions  (and  the  initial  condition  appears  in  the  last  equation  for  e). 
Once  these  boundary  conditions  are  specified,  the  components  of  6  are 
determined  by  decomposition  to  some  n-term  approximation.  These  conditions 
are: 

0X  =  0  at  x  =  0  and  x  *  % 

0 ( t  =  0)  is  a  specified  function  f(x,y) 

9y(y  =0)  is  a  specified  function  g(x) 

l9y  +  B9la  =  0 

where  a  is  a  constant. 

Decomposition  calculations  are  simpler  than  by  other  methods,  and,  most 
important,  solutions  especially  in  nonlinear  or  stochastic  cases  become 


more  realistic.  Of  course,  modelling  is  necessary  before  we  can  apply 
any  method  and  we  must  begin  with  a  model.  Now  that  we  know  decomposition 
works,  we  need  not  begin  with  simplistic  models. 


14)  Transient  heat  ■ conduction :  An  equation  describing  transient  heat 
conduction  is  a  region  R  is 

V^u  +  k”V(x,y,z,t)  a  a-13u/3t 

u  is  temperature,  k  is  thermal  conductivity,  and  a  is  thermal 
diffusivity.  k  and  a  are  constants.  Six  boundary  conditions  and  an 
initial  condition  are  required.  (Such  a  problem  is  described  by  Beck  in 
AIAA  Journal ,  24,  no.  2,  Feb.  1984,  pp.  327-333  with  the  use  of  Butkovsky's 
Green's  functions?)  In  our  form,  we  write 


CL*x  +  Lyy  +  4*1“  +  k',f  "  “'V 


We  solve  solve  for  the  four  linear  operator  terms,  apply  the  appropriate 

00  oo 

inversions,  let  u  =  Y  u  and  f  =  Y  A„  where  the  A„  are 

n=0  n  n=0  n 

generated  for  the  function  f. 


L„„u  =  a 

’^L.u 

-  L 

u 

-  L  u  -  k”^ 

XX 

t 

yy 

zz 

1 _ u  =  a 

-  L 

u 

-  L  u  -  k"1 

yy 

t 

ZZ 

XX 

L  u  =  a 

-\u 

-  L 

u 

-  L  u-k1 

zz 

t 

XX 

yy 

L.u  =  aL 

u  + 

aL 

u 

+  aL  u  +  k"^ 

t 

XX 

yy 

zz 

A.  G.  Butkovsky,  Green's  Functions  and  Transfer  Functions  Handbook, 
Wiley,  1982. 
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u  =  d>  +  a-1L_1Uu  -  L-1L  u  -  L-1L  u  -  k'W 
*xx  xx  t  xx  yy  xx  zz  xx 

U  =  d>  +  a"^L_1L.u  -  L_1L  u  -  L_1L  u  -  k"Vlf 
vyy  yy  t  yy  ZZU  yyuxx  yy 

u  -  <|>zZ  +  a  LzzLtu  “  LzzLxxu  ’  LzzLyyu  "  k  Lzzf 
u  =  d>t  +  aL;]Lxxu  +  aLj\u  +  “L‘1Lzzu  +  k^L^f 


‘t  ‘-yy  t  uzzv 


<j>t  satisfies  the  initial  condition  u(x,y,z,0)  and  the  4>xx ,  <j>^,  <t>zz 
each  satisfy  two  boundary  conditions.  The  <f>xX,  <j>  ,  <J>ZZ,  4>t  are  the  uQ 

oo 

terms  of  the  decomposition  for  each  equation.  Then  since  u  =  T  u 

n=0  n 

OO 

and  f  =  l  A  the  following  components  are  determined.  For  example 
n=0  n 

U1  "  “  LxxLtut  "  LxxLyyu0  "  LxxLzzu0  "  k  LxxA0 

U1  ”  “  LyyLtu0  "  LyyLzzu0  "  LyyLxxu0  "  k  LyyA0 

U1  =  °  ^LzzLtu0  “  LzzLxxu0  *  LzzLyyu0  "  k  1lzzA0 

U1  =  aLt^Lxxu0  +  alVLyyu0  +  aLtlLzzu0  +  k  W 

etc.,  tor  U2»  u^,  ...  .  The  solution  u  is  obtained  by  adding  the  four 
equations  for  u  and  dividing  by  4,  or,  as  an  approximation,  by  using 
the  four  approximate  solutions  4>n . 
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NONLINEAR  OSCILLATIONS 


Introduction:  Practically  all  the  problems  of  mechanics  are 

nonlinear  at  the  outset.  Nonlinear  oscillating  systems  are  generally 
analyzed  by  approximation  methods  which  involve  some  sort  of  linearization. 
These  replace  an  actual  nonlinear  system  with  a  so-called  "equivalent" 
linear  system  and  employ  averagings  which  are  not  generally  valid. 

While  the  linearizations  commonly  used  are  adequate  in  some  cases, 
they  may  be  grossly  inadequate  in  others.  Thus,  if  we  want  to  know  how 
a  physical  system  behaves,  it  is  essential  to  retain  the  nonlinearity, 
not  just  solve  a  convenient  mathematized  model. 

Using  the  decomposition  method,  restrictive  assumptions  on  nonlinear 
and  stochastic  behavior  for  the  sake  of  a  tractable  model,  become 
unnecessary;  correct  solutions  are  obtained  in  the  strongly  nonlinear 
case  and  in  the  case  of  stochastic  (large  fluctuation)  behavior,  as  well 
as  in  the  cases  where  perturbation  would  be  applicable  or  in  the  linear 
and/or  deterministic  limits. 

We  are  concerned  in  this  section  with  the  study  of  vibrations  or 
oscillatory  motion  and  the  associated  forces.  Vibrations  can  occur  in 
any  mechanical  system  having  mass  and  elasticity.  Consequently,  they 
can  occur  in  structures  and  machines  of  all  kinds.  In  large  space 
structures  containing  men  or  machinery,  such  vibrations  will  result 
in  difficult  and  crucial  control  problems  and  also  lifetime  or  duration 
considerations  since  vibrations  can  lead  to  eventual  failure. 

Oscillations  can  be  regular  and  periodic,  or  they  can  be  random. 
Randomness  leads  to  stochastic  differential  equations.  In  deterministic 


systems  -  the  special  case  where  randomness  vanishes  -  the  equations  modeling 
the  system  provide  instantaneous  values  for  any  time.  When  random  functions 
are  involved,  the  instantaneous  values  are  unpredictable,  and  it  is 
necessary  to  resort  to  a  statistical  description.  Such  random  functions 
of  time,  or  stochastic  processes,  occur  in  problems,  for  example,  such  as 
pressure  gusts  encountered  by  aircraft,  jet  engine  noise,  or  ground  motion 
in  earthquakes. 

Oscillatory  motion  is  modelled  by  equations  of  the  general  form: 

y"  +  f(y.y’.t)  =  x(t) 

Fy  =  x ( t) 

Stochastic  processes  may  be  involved  in  coefficients,  input  terms,  or 
initial  boundary  conditions,  so  x(t)  can  be  assumed  to  be  generally 
stochastic.  Input  conditions  -  possibly  stochastic  -  are  given  (statistically 
described  if  stochastic).  Therefore  F  will  be  a  nonlinear  stochastic 
operator  with  (possibly)  stochastic  coefficients.  More  conveniently,  we 
write  Fy  =  Ly  +  Ny  where  Ly  is  a  linear  (stochastic)  term,  and  Ny 
is  a  nonlinear  (stochastic)  term,  or 

Ly  +  Ny  =  x. 

Now,  for  the  linear  operator  i,  define  L  =  L  +  R  +  R  where:  L  is 
a  linear  invertible  deterministic  operator,  so  that  L”1  exists;  R  is 
the  remainder  of  the  linear  deterministic  operator,  vanishing,  of  course, 
if  we  can  easily  invert  the  entire  linear  deterministic  operator,  and  R 
is  a  (linear)  stochastic  operator. 
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Similarly  for  the  nonlinear  operator  N ,  define  Ny  =  Ny  +  My  since 
there  may  be  both  deterministic  nonlinear  terms  Ny  and  stochastic 
(nonlinear)  terms  My. 

Thus  in  general  we  rewrite  our  general  oscillator  equation  as 


Ly  +  Ry  +  i?y  +  Ny  +  wy  =  x. 


Note  that  in  a  particular  problem  any  number  of  terms  from  one  to  four  may 

vanish.  Still  more  generally,  Ny  may  actually  be  a  function  of  y.y',... 

but  we  will  again  write  Ny  for  simplicity.  Calculation  of  the  An  for 

such  a  case  becomes  more  difficult  and  will  not  be  discussed  here.  Our 

purpose  now  is  not  discussion  of  the  entire  mathematical  theory  but 

demonstration  of  applicability  to  space  structures. 

Convergence  of  the  decomposition  series  can  reasonably  be  expected  to 

be  fastest  when  we  invert  the  entire  linear  deterministic  operator. 

Computation  of  the  integrals  will,  however,  be  more  difficult  should  this 

be  the  case,  since  the  Green's  functions  will  not  then  be  simple.  Thus, 

for  ease  of  computation  we  will  let  L  denote  the  highest  order  differential 
2  2 

operator  or  d  /dt  in  the  above  oscillator  equation. 

Let  us  review  the  significance  of  the  mathematical  expressions  we 
will  deal  with  in  order  to  make  the  work  as  clear  as  possible.  In  an 
oscillator,  we  have  generally  an  external  force  or  driving  term  x(t),  a 
restoring  force  f(y)  dependent  on  the  displacement  y,  and  a  damping 
force,  since  energy  is  always  dissipated  in  friction  or  resistance  to 
motion.  Usually  this  is  dependent  on  velocity,  and  we  will  write  it  as 


g(y'). 


If  we  have  a  free  oscillating  mass  m  on  a  spring  with  no  damping, 


we  can  write 


WDTOT".V'Jr». 


my"  +  ky  =  0 


if  the  spring  obeys  Hooke's  Law,  i.e.,  assuming  displacement  proportional 
to  force.  Of  course,  no  spring  really  behaves  this  way.  Often  the  force 
needed  for  a  given  compression  is  not  the  same  as  for  an  extension  of  the 
same  amount.  Such  asymmetry  is  represented  by  a  quadratic  force,  or  force 

p 

proportional  to  y  rather  than  y.  We  may  have  a  symmetric  behavior 
with  proportionality  to  y  .  In  this  case,  while  the  solution  is  not  the 
harmonic  solution  which  one  gets  for  the  model  equation  my"  +  ky  =  0, 
it  is  still  a  periodic  solution.  The  damping  force  g(y')  may  be  cy' 
where  c  is  constant,  or  it  may  be  more  complicated  such  as  a  nonlinear 
function  of  y  or  its  derivatives,  e.g.,  g(y',y'  )  so  that  it  depends 
on  v  as  well  as  v.  By  usual  methods,  analytic  solutions  then  become 
impossible,  i.e.,  we  can  not  deal  analytically  with  such  a  case  unless  we 
resort  to  linearization. 

If  we  write  -f (y )  for  the  restoring  force,  -h(y')  for  the  damping 
force,  and  represent  the  driving  force  with  g;  the  resulting  equation 


will  be 


y"  +  f(y)  +  h(y')  =  g 


Suppose  the  restoring  force  is  represented  by  an  odd  function  such  that 
f(y)  =  -f(-y).  We  have  this  in  most  applications;  it  means  simply  that 
if  we  reverse  the  displacement,  then  the  restoring  force  reverses  its 
direction.  A  pendulum,  for  example,  behaves  this  way.  We  might  take  the 
first  two  terms  of  the  power  series  for  f(y)  and  write  f(y)  =  ay  +  By  ■ 


Then  we  have 


y"  +  ay  +  By3  =  g 
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If  we  have  damping  also,  we  have 

3 

y"  +  cy'  +  ay  +  0y  =  g 

assuming  the  damping  force  is  -cy1.  (This  is  Duffing' 8  equation.) 

We  will  begin  -  after  a  simple  example  of  a  pendulum  (harmonic 
motion)  -  with  the  well-known  anharmonic  oscillator  then  go  on  to 
consider  more  general  oscillators  such  as  the  Duffing  oscillator  and 
the  Van  dev  Pol  oscillator.  The  Duffing  oscillator  in  a  random  force 

O 

field  is  modeled  by  y"  +  ay*  +  By  +  y y  =  x(t).  It  can  be  analyzed 
without  limiting  the  force  x(t)  to  a  white  noise  or  restricting  a,  6, 
y  to  be  deterministic.  The  same  applies  to  the  Van  der  Pol  oscillator 

o 

modeled  by  y"  +  ey y'  -  ey'  +  y  =  x(t).  These  equations  are  in  our 
standard  from  Fy  =  x(t)  which  can  be  solved  by  the  decomposition 
method.  If  the  equation  is  linear  and  deterministic,  we  have  Fy  =  Ly  = 
An  equation  that  is  deterministic  but  nonlinear  can  be  rewritten  as 
Fy  =  Ly  +  Ny  =  x.  A  linear  stochastic  equation  is  Ly  =  x,  etc.  These 
cases  are  discussed  in  the  earlier  work.  Consider  the  pendulum  problem 
as  an  example,  then  we  proceed  with  the  anharmonic  oscillator. 

A  Simple  Nonlinear  Oscillator  Example: 

Consider  the  simple  vertical  pendulum  consisting  of  a  mass  m  at 
the  end  of  a  rod  length  l  moved  through  an  angle  0.  We  obtain 
immediately 

d20/dt2  +  k2sin  0=0 

where  k2  =  g/2..  Let  L  =  d2/dt2  and  N(0)  =  k2  sin  B  to  obtain  our 
usual  standard  form  L0  +  N0  =  0  for  a  homogeneous  deterministic 
differential  equation. 
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Small-amplitude  case:  The  usual  treatment  is  to  simplify  this 

2  2  2 

anharmonic  motion  by  assuming  sin  0  *  0.  We  can  then  write  d  0/dt  +  k  0  =  0, 
the  well-known  harmonic  oscillator  problem.  L  is  defined  as  before,  but 
now  there  is  no  nonlinear  term,  and  the  remainder  of  the  linear  term  is  R0, 
i.e.,  the  entire  linear  operator  is  decomposed  into  L  +  R  where  R  is 
a  linear  deterministic  "remainder"  operator.  Define  L“^  as  the  double 
definite  integration  from  0  to  t.  If  we  have  a  forcing  term  as  well, 
we  would  have  L0  +  R0  =  x(t).  The  solution  is  found  by  writing 


L0  =  x  -  R0 


L_1L0  =  0  -  0(0)  -  t0 1 (0)  =  L'^x  -  L_1R0 


0(t)  =  0Q  -  l  R0 

with  0O  =  0(0)  +  t0'(O)  +  L_1x(t), 


Now  substituting  0  =  l  en(t). 


we  find 


Vi  *  -L'lRen 


so  that  all  components  can  now  be  determined. 

Let  us  assume  initial  conditions  are  given  as  0(0)  =  y  and 
6'(0)  *  0.  Then,  since  x  =  0 


e0  =  Y 

e1  =  -L_1R0o  =  -L'Vy  =  -ykZtZ/2! 

02  =  -L"1R0]  =  -L-1k2(-Yk2t2/2! )  =  yk4t4/4I 

03  =  -yk6t6/6: 
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We  have 


6n  =  Y(-l)n'1(kt)2n_2/(2n-2): 

6=y  l  (-l)n_1(kt)2n'2/(2n-2): 
n=l 


Thus  6  =  y  cos  kt  with  k  =  (g)  '  . 

If  we  assume  initial  conditions  8(0)  =  0,  6 ' (0)  =  w,  we  obtain 


0q  =  tat  -  (u)/k)kt 
e]  =  -(w/k)k3t3/3! 
02  =  (oj/k)k5t5/5! 


6n  =  (w/k)  (-1  )n_1  (kt)2n~V(2n-l )! 


and  finally. 


e(t)  =  (w/k ) tkt  -  (kt)3/3:  ♦  (kt)5/s:  - 
■  (Wk)  I".,  (-1  )n_1  (kt)zn_1/(2n-1 ): 


0(t)  =  (w/k)  sin  kt 


The  case  6(0)  =  y»  and  8'(0)  =  w  yields 


6q  =  y  +  (w/k)kt 

6-j  =  -L_1R[y  +  (u/k)kt] 


e„  =  Y(-l)n_1  (kt)2n_2 /  ( 2n-2 ) - 


+  (w/k)(-l)n_1  (kt)2n_1/(2n-l): 


ixvwvTrv  iru  jtujtvj 
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or 

0(t)  =  y  cos  kt  +  (w/ kj  sin  kt 
the  well-known  general  solution  of  the  harmonic  oscillator. 

Large -amplitude  case:  We  can  generalize  to  the  anharmonic  oscillator 
describing  the  pendulum  problem  for  large  amplitude  motion,  the  only 
difference  being  that  since  we  now  have  a  nonlinearity,  it  must  be 
expressed  in  terms  of  the  polynomials. 

Consider  the  equation  now  for  large-amplitude  motion  in  the 
pendulum 

d2e/dt2  +  k2  sin  9  =  0 

with  k2  *  g/Jl.  Let  L6  *  d2e/dt2  and  N(0)  *  k2  sin  0.  We  let  the 
nonlinear  term  N0  be  represented  by  the  sum  of  the  An  polynomials 
generated  fro  the  nonlinear  term  as  usual  in  the  decomposition  method. 


Thus  N0  =  7  A„  where  the  A„  are  given  by: 
n=0  n  n 


A0  =  sin  e0 


A-j  =  01  cos  0q 


A2  =  -(0^/2)  sin  9q  +  62  cos  0O 


A3  =  -(0^/6)  cos  0Q  -  0-|02  sin  0Q  +  83  cos  0Q 
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This  equation  is  a  nonlinear  deterministic  homogeneous  second-order 
differential  equation.  We  must  assume  appropriate  initial  conditions. 
We  will  choose  here  8(0)  =  y  =  constant  and  8'(0)  =  0.  We  can  do  it 
just  as  well  for  other  combinations  of  initial  conditions,  of  course, 
as  will  be  obvious.  We  now  have 


16  =  -N(9) 


Operating  with  the  inverse  of  L,  a  two-fold  integration,  we  have 
0  *  0O  -  L_1N(0) 

where 

0Q  =  6(C)  +  te'(o) 

which  in  our  case  is  just  y.  The  N(8)  term  -  we  will  henceforth 

00 

drop  the  parentheses  -  is  replaced  by  the  J  A  .  We  now  have 

n=0  n 

00 

e'Y-L'1  l  An(k2  sin  0) 
n=0  n 

We  prefer  to  write 

0  =  y  -  L'V  l  A  (sin  6) 
n=0  n 

2  2 

i.e.,  we  can  let  N0  =  sin  0  rather  than  k  sin  0  since  k  is  a 
constant.  The  corresponding  An(sin  0),  or  simply  An>  are  [1]: 


A0  =  sin  e0 
A^  =  0^  cos  0Q 


A2  =  -(02/2)  sin  0Q  +  02  cos  0Q 

A3  =  - (0 -j /6 )  cos  0Q  -  8^2  sin  0Q  +  03  cos  0Q 
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etc. ,  to  yield 


•t  ft 
).  =  -  dt 
1  ■’0  J0 


dt  k2  sin  0q  = 


t  ft 


fi 

dt 

JO  J0 


2  . 

dt  k  sin  y 


=  -(sin  Y)(k2t2/2!) 


=  -  |  dt  |  dt  k2(e1  cos  eQ) 


’0  ■'0 


ft  ft 

dt  dt  k2[-(k2t2/2!)sin  y]  cos  y 
JO  h 


=  ( k^t^/4 ! ) sin  y  cos  y 


e3  =  -  |  dt  j  dt  k2[-(e2/2)sin  eQ  +  e2  cos  eQ] 


'0  '0 


-(k6t6/6! ) * [si n  y  cos2  y  -  3  sin3  y] 


etc.  Thus, 


e(t)  =  y  -  [(kt)2/21]sin  y  +  [(kt)4/4!]sin  y  cos  y 


-  [(ktr/6!  )][sin  y  cosy  -  3  sin^  y] 

+  [(kt)®/8!][-33  sin3  y  cos  y  +  sin  y  cos3  y] 


l.e.,  we  have  the  solution  for  large-amplitude  motion. 
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As  a  check  we  can  let  y  be  sufficiently  small  so  that  small  amplitude 


motion  is  being  considered.  Then 

e(t)  =  y[1  -  (kt)2/2!)  +  (kt)4/4!  -  •••] 

which  is,  of  course,  the  result  for  the  linear  harmonic  oscillator  with 
the  given  initial  condition. 


IS)  The  Duffing  and  Van  dev  Tot  Oscillators:  In  the  preceding  section  the 

decomposition  method  was  applied  to  harmonic  and  anharmonic  oscillators. 

We  now  apply  it  to  the  Duffing  oscillator  and  the  Van  der  Pol  oscillator. 

It  requires  no  linearization  or  "smallness"  assumptions.  The  treatment 

can  also  include  randomness  in  coefficients  or  inputs  without  customary 

restrictions  to  special  processes  or  perturbation  theory. 

When  we  consider,  for  example,  the  equation  for  a  simple  pendulum,  we 

approximate  sin  x  by  x  to  obtain  the  harmonic  oscillator  equation. 

3 

Suppose  we  go  a  step  further  and  write  sin  x  =  x  -  x  /3!,  i.e.,  use  the 

first  two  terms  of  the  series  for  sin  x  assuming  small  x.  We  then 
get  the  equation  x"  +  w  x  +  ex  =  0,  which  is  the  Duffing  equation  with  e 
as  a  "small"  parameter.  This  is  an  example  of  a  perturbation  method;  one 
would  seek  a  solution  in  the  form  x(t)  =  xQ(t)  +  ex^(t)  +  ...  .  Using 
decomposition,  we  are  not  restricted  to  assuming  small  parameters. 

Statistical  Linearization:  In  dealing  with  stochastic  oscillators,  we 
depart  again  from  usual  procedures  which  require  some  sort  of  approximation 
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in  order  to  determine  the  second-order  response  statistics.  A  common 
procedure  in  this  connection  is  statistical  linearization.  This  procedure 
simply  replaces  the  original  nonlinear  equation  with  a  so-called  "equivalent" 
linear  system.  Thus,  if  we  write  an  oscillator  equation  in  the  form: 

x"  +  ax'  +  UqX  +  8f(x)  =  F(t) 

where  x(t)  is  a  displacement,  a  is  a  damping  constant,  Wq  is  a 
linear  frequency,  8f(x)  is  a  nonlinear  restoring  force,  and  F(t)  is 
a  stationary  process,  the  process  of  statistical  linearization  substitutes 

x"  =  ax'  +  y^x  =  F(t) 

2 

where  y  is  determined  in  such  a  way  that  the  mean  square  error 
due  to  the  replacement  is  minimized,  and  the  mean  displacement  is  the 
same  for  both  systems.  It  is  customary  to  assume  F(t)  is  Gaussian  and 
delta-correlated  with  zero-mean,  or,  <F(t)>  =  0  and  <F(t)F(t')>  -  DS(t-t'). 
This  latter  assumption  is,  of  course,  made  for  mathematical,  not  physical, 
reasons  and  is  physically  unrealistic.  We  propose  none  of  these 
restrictions  and  will  solve  the  actual  nonlinear  equation. 


The  Decomposition  Method:  The  Duffing  Oscillator  is  described  by 
the  equation: 

y"  +  ay’  +  By  +  yy3  =  x(t) 

In  our  standard  Fy  =  x(t)  form. 


The  Van  der  Pol  equation  is  generally  given  as: 


or  by 


y"  +  ey2y'  -  ey'  +  y  *  x(t) 


y"  +  ey'(y2  -  1)  +  y  =  x(t) 


which  we  can  write  in  the  form 


y"  +  ay'  +  By  +  y(d/dt)yJ  =  x(t) 

since  y2y'  =  (d/dt) (y^/3) .  Thus  a  =  -e,  B  =  1.  Y  =  e/3  relates 
the  last  equation  to  the  two  previously  given  forms.  We  now  have  our 
standard  form  Fy  =  Ly  +  Ny  =  x,  or  Fy  =  Ly  +  Ny  =  x  (if  no  stochastici ty 
is  involved).  Here  L  is  assumed  to  be  a  linear  and  invertible  operator 
and  N  is  a  nonlinear  operator.  We  will  consider  the  equations  to  be 
deterministic  here. 

The  linear  operator  L  in  both  the  Duffing  and  Van  der  Pol  oscillator 

equations  is  given  by  d/dt  +  ad/dt  +  B.  The  nonlinear  term  Ny  is  a 

3 

simple  cubic  nonlinearity  yy  in  the  case  of  the  Duffing  oscillator,  and 
y(d/dt)y  in  the  case  of  the  Van  der  Pol  oscillator.  These  terms  will, 
as  usual,  be  expanded  in  our  An  polynomials  generated  for  the  specific 
nonlinearity. 

The  treatment  of  the  linear  operator  offers  some  alternatives.  We 
can  use  the  entire  linear  operator  as  L  which  enhances  speed  of 
convergence,  but  the  inverse  and  consequent  integrations  become  more 
difficult.  We  can  also  use  part  of  the  above  operator  which  could  be 
L  -  d2/dt2,  L  =  d2/dt2  +  ad/dt,  or  L  =  d2/dt2  +  B.  We  prefer  in  most 

O  p 

cases  to  use  L  =  d  /dt  ,  i.e.,  the  highest  order  differential  operator. 
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We  expect  this  to  give  the  slowest  convergence  but  much  easier  integrations 
and  less  actual  computation  time. 

The  remainder  of  the  linear  operator  will  be  called  R,  the  "remainder" 

2  2 

operator.  If  L  =  d  /dt  ,  R  =  ad/dt  +6.  (If  we  also  consider  stochasticity, 
we  will  use  a  script  letter  r  for  a  random  part  of  the  operator  and  may 
have  L  +  R  +  R.) 

The  choice  made  here  (that  L  =  d  /dt  )  yields  the  simplest  Green's 
function  for  computation.  In  this  case,  L"^  is  the  two-fold  definite 
integral  from  0  to  t.  Generally,  this  choice  of  the  highest  ordered 
derivative  for  L  is  the  most  desirable  because  the  integrations  are 
the  simplest.  If  we  invert  the  entire  linear  operator,  convergence  is 
expected  to  be  much  faster.  It  is  interesting  to  examine  a  compromise 
here  which  can  be  used  to  advantage  on  occasion. 

If  we  choose  L  =  d2/dt2  +  B,  R  a  ad/dt,  we  gain  something  in 
convergence  rate  over  the  previous  case  and  expect  to  lose  something 
in  ease  of  computability.  The  interesting  aspect  that  suggests  the 
compromise  is  that  we  see  we  will  get  sine  and  cosine  functions  for 
solutions  of  the  homogeneous  equation.  For  the  Duffing  equation  we  now 
have 

Ly  =  x  -  Ry  -  yy3 

y  =  C^U)  +  c2<J>2(t)  +  L_1x  -  L"1  Ry  -  YL_1y3 
where  ^ ,  <f>2  satisfy  L<j>  =  0  or  d^/dt2  +  6<f>  =  0.  Consequently, 

4>}(t)  =  cos/Bt 


^(t)  =  (l//f?)  sin  /Bt 
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w  oo 

l  yn  =  c14>1(t)  +  Cp^(t)  +  L"x  -  L-1Ry  -  yl~]  l  An 
n=0  "  c  c  n=0  n 


I  yn  =  c,  cos  /gt  +  (c9//g)  sin  /gt  +  L-1x 
n=0  n  2 

-  aL_1(d/dt)  l  y  -  yL_1  £  A 
n=0  n  n=0  n 

where  the  An  are  the  appropriate  polynomials  for  Ny  =  y3.  These 
are  given  by:  -  . 

A0  =  y0 

A1  =  ^1 

A2  ■  3y„y3  +  3y^y2 

V'l*  6Wly2  *  3y0y3 

\  =  3yly2  +  3yOy2  +  6y0yly3  +  3y0y4 

As  =  3y,yf  +  3y3y3  ♦  6yoy2y3  *  6y0y,y4  +  3y3y5 

A6  *  y|  +  6yly2y3  *  3yly4  +  3V3  +  6y0y2y4  +  6Vly5  +  3y0y6 


2  2  -1 

Since  L  =  d  /dt  +  B  now,  L  is  no  longer  the  simple  two-fold 
integral,  and  we  must  determine  the  Green's  function  for  this  L.  The 
G  will  satisfy  the  equation  LG(t,x)  =  <5(t  -  x)  or 

d2G(t,x)/dt2  +  0G(t,x)  =  6(t  -  x) 
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G,  of  course,  is  determinable  in  a  number  of  ways.  We  will  again  use 
the  decomposition  method  itself  and  so  we  prepare  by  writing 


d2G/dt2  =  6 ( t  -  t)  -  $G 


Again  we  have  a  simple  second-order  operator  to  invert.  Hence 

00 

G(t.-r)  =  G(0,t)  +  tG.  (0,t)  +  L_16(t  -  t)  -  6L_1  l  G 

1  r>  =  n  " 


OO  OO 

J.  Gn  =  G(0,t)  +  LG.(O.t)  +  L-16(t  -  t)  -  8L'1  l  G 
n=0  n  1  n=0  n 


■  G0  -  eL  l  Gn 
u  n=0  n 


where 


Gq  =  G(0,t)  +  Gt(0,t)  +  tH(t  -  t) 

G1  -  -3L‘]G0  =  -p{G(0,T)t2/2:  +  Gt(0,T)t3/3! 


+  (t3/3! ) H ( t  -  t)> 


•BG(0,T)[t2/2:  +  t3/3!]  -  B(t3/3:)H(t  -  t) 


etc.,  for  G3,  ...  .  An  appropriate  n-term  approximation  can  now  be 

used  in  the  l”1  integrations. 

2  3 

The  example  my"  +  n  y  +  ay  =0  occurs  in  the  theory  of  nonlinear 
vibrating  mechanical  systems  and  in  some  nonlinear  electrical  systems 
which  is  the  Duffing  case  above  (with  no  damping).  Suppose  we  have  the 
specified  conditions  y  =  a  at  t  =  0  and  y'  =  0  at  t  =  0.  Write 
L  =  d2/dt2  and 
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Ly  =  -(w2/m)y  -  (a/m)y3 


y  =  y(0)  -  L"1  (u)2/m)  l  y  -  (a/m)  £  A 

n=0  n  n=0  1 

where 

*0  =  3 

y1  =  -(u2/m)L_1y0  -  (o/m)A0 

2  2  3 

=  -aw  t  /2m  -  aa  /m 

so  that  y  -  a[l  -  w2/t2/2m  -  aa2/m  . ..]. 


Solution  of  the  General  Case  of  Duffing’ s  Equation:  For  the  Duffing 
2  2  °° 

equation  with  L  =  cr/dt^  we  have  y  =  £  y  where 

n=0  n 

y0  =  y(0)  +  ty'(0)  +  i-'Mt) 

y-j  =  -L_1a(d/dt)y0  -  L'Vq  -  L'Vg 

y2  =  -L"1a(d/dt)y1  -  L'1By1  -  L"1yA] 

y3  «  -L_1a(d/dt)y2  -  L_18y2  -  L_1yA2 


Solution  of  Van  der  Pol  equation:  For  the  Van  der  Pol  equation,  we 


have 

y0  =  y( o)  +  ty'(o)  +  L-1x(t) 

y]  =  -L_1a(d/dt)y0  -  L_1By0  - 

y2  =  -L'1a(d/dt)y1  -  L"16y1  -  L'1y(d/dt)A1 

etc. 


16)  Stochastic  Oscillations:  We  could  also  have  stochastic  fluctuations 
in  a,  8,  or  y  in  addition,  of  course,  to  stochastic  x(t)  or  initial 
conditions.  Thus,  in  general  we  could  write 

a  =  <a>  +  e 

6  =  <8>  +  n 
y  *  <y>  +  o 

where  e,  n,  o  are  zero-mean  random  processes.  The  solution  process  can 
now  be  obtained  from 

oo 

Ly  =  x  -  a(d/dt)y  -  By  -  y  l  A_ 

n=0  n 

OO 

-  e(d/dt)y  -  ny  -  a  7  A„ 

n=0  n 
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3  3 

where  the  An  summation  represents  y  in  the  Duffing  case  and  (d/dt)y 

in  the  Van  der  Pol  case.  Thus, 

y0  =  y(0)  +  ty'(O)  +  L-1x 
y]  =  -a(d/dt)yQ  -  ByQ  -  yAq 

-  e{d/dt)yQ  -  nyQ  -  oAQ 

y2  =  -a(d/dt)y-|  -  -  yA] 

-  e(d/dt)y1  -  ny1  -  oA-j 

Thus  the  summation  of  the  above  terms  to  some  yp  is  the  solution  of  both 
oscillator  equations  except  the  A  are  different,  of  course,  as  explained 

oo 

above.  Then  y(t)  =  £  y  (t)  yields  a  stochastic  series  from  which 

n=0  n 

statistics  can  now  be  obtained  easily  by  averaging  without  problems  of 
statistical  separability  of  quantities  such  as  <Ey>  where  E  =  e(d/dt)  -  n 
which  normally  require  closure  approximations  and  truncations. 

Eernark:  When  transient  behavior  becomes  oscillatory  and  periodic 
after  a  certain  time,  we  can  then  consider  it  as  a  boundary-value  problem. 
Sometimes  one  can  obtain  a  solution  for  the  oscillating  behavior 
directly  starting  with  the  transient  behavior.  In  other  cases  it  may 
have  to  be  approached  separately. 


COMPUTATIONAL  ASP:CTS  OF  STRUCTURAL  MECHANICS 


These  have  been  recently  discussed  by  Noor  and  Atluri.  In  proposed 
large  flexible  space  structures,  supercomputer  simulations  will  be 
important  but  are  likely  to  be  so  large  and  complex  that  any  means  of 
reducing  the  computational  effort  needs  to  be  a  paramount  objective.  The 
discretization  inherent  in  the  computer  techniques  makes  it  difficult  to 
see  dependences  or  gain  insight  into  the  nature  of  the  response  and, 
further,  enormously  increases  computation  as  grid  spacing  is  decreased 
for  accuracy. 

The  decomposition  method  provides  continuous  rather  than  discretized 
solution  -  a  solution  in  analytic  form  -  with  an  enormous  potential 
saving  in  computation.  This  feature,  as  well  as  the  natural,  and  physical, 
treatment  of  nonlinearity  ar.d  stochasticity  in  the  decomposition  method, 
and  its  computability  provide  a  powerful  new  tool  for  structural  analyses. 


17)  Stochastic  Structural  Dynamics:  The  equation  Fu  -  g(t)  can  represent 
a  generic  problem  in  vibration  where  F  is  a  nonlinear  stochastic 
operator  and  g  is  a  stochastic  process.  Suppose  Fu  =  Lu  +  tfu  where  L 
is  a  linear  stochastic  operator  and  nu  is  a  nonlinear  (and  possibly 
stochastic)  term.  L  can  be  further  decomposed  into  L  +  R  +  R  where 
L  is  a  deterministic  operator-specifically  the  highest-order  linear 
differential  operator  -  typically  d  /dt  in  a  vibration  problem,  R  is 
the  "remainder  operator"  representing  the  remaining  portion  of  the  linear 
and  deterministic  operator.  (This  device  removes  the  necessity  of 


i'  i.lt.iMUt.n 


determination  of  the  Green's  function  for  the  entire  linear  deterministic 


operator),  and  r  indicates  a  stochastic  operator. 


If  we  consider,  for  example,  the  equation 


my( t)  +  a  y(t)  +  b  y(t)  =  f(t) 


and  divide  through  by  m,  we  can  write  now 


y(t)  +  a  y(t)  +  0  y(t)  =  g(t) 


2  2 

Let  L  =  d  /dt  ,  R  =  a  d/dt.  Now  we  will  assume  for  an  example  that  a 


is  constant  but  3  is  a  random  process  B(t)  which  quite  reasonably,  is 


assumed  statistically  independent  of  the  input  process  g(t).  In  the 


decomposition  notation  consistent  with  the  notation  in  the  literature, 


we  let  R  -  B(t).  We  have 


Ly  +  Ry  +  Ry  =  g(t) 


Ly  =  g  -  Ry  -  Ry 


yn  =  A  +  Bt  +  L"  g 


-L'V  '  L'Vn 


yg  is  known  as  soon  as  specific  initial-boundary  conditions  are  given. 
All  following  components  are  determined  for  n  >  0.  Now  y  =  yn 


and  an  n-term  approximation  is  given  by  <fc  =  y^ .  Now  the 
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expectation  of  y  Is  approximated  by  <$n>.  Second-order  statistics 
are  also  determinable  -  the  correlation  KU-j.tg)  is  given  by 
<4>n(ti )<j>n(t2)>-  This  method  avoids  all  the  usual  restrictive  assumptions 
and  there  is  no  statistical  separability  problem  which  ordinarily  leads 
to  closure  approximations  in  averaging  or  hierarchy  methods.  Conventionally, 
stochastic  structural  dynamics  can  only  allow  the  system  input  g  to  be 
stochastic.  The  decomposition  method  treats  linear  and  nonlinear  equations 
without  closure  approximations  or  perturbation.  In  addition,  randomness 
can  be  considered  in  parameters  or  conditions.  These  cases  are  simply 
described  in  [2]. 

Orada  and  Haftka  point  out  (in  the  AIAA  Journal,  vol.  25,  no.  8, 

Aug.  1987,  pp.  1133-1138)  that  in  active  vibration  control,  because  of 
the  interaction  between  the  structure  and  the  control  system,  simultaneous 
optimal  design  of  both  systems  may  be  necessary.  (More  traditionally, 
the  structure  is  optimized  to  minimize  weight  subject  to  the  stress  and 
stiffness  constraints  while  the  control  system  is  optimized  to  minimize  a 
performance  measure  involving  deformation  and  control  effort.)  Such 
simultaneous  optimization  is  discussed  in  the  referenced  paper  using  white 
noise  inputs  and  linear  control  theory.  We  would  point  out  only  that  a 
more  powerful  tool  has  now  become  available  for  such  optimizations.  We 
have  seen  now  that  the  control  theory  can  be  made  nonlinear  and  stochastic 
through  application  of  the  decomposition  method  without  the  use  of 
numerical  methods  and  we  will  be  able  to  pursue  this  when  funding  is 
available. 

The  decomposition  method  allows  more  realistic  models,  eliminating 
restrictive  mathematical  assumptions.  As  an  example,  consider  the 
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equation:  x  +  ax  +  0x  =  g  where  vibrating  motion  in  the  x  direction,  ax 

is  the  damping  force  and  0x  is  the  restoring  force.  Initially  the  system  is 
at  rest  so  x(0)  *  x(0)  =  0.  We  can  assume  a,  0,  g,  or  any  combination,  to 
be  stochastic.  The  widely-used  averaging  or  hierarchy  methods  of  determining 
statistics  (seeking  moments  of  the  solution  process  directly)  lead  to 
unjustified  closure  approximations.  Using  decomposition,  we  do  not  need  to 
require  stationarity  which  is  physically  unlikely  and  mathematically  restrictive. 
Perturbation  methods  work  best  in  cases  where  fluctuations  are  small.  In  the 
perturbative  case,  one  is  assuming  corrections  to  a  deterministic  solution  are 
minor. 

In  a  space  borne  effort  costing  many  billions  of  dollars,  it  is 
essential  that  we  seek  physical  solutions  rather  than  mathematical  solutions 
to  models  assuming  processes  that  either  do  not  exist  in  nature  or  exist 
in  cases  of  minor  vibrations.  The  mathematics  must  take  into  account  any 
possibility  of  vibrations  or  effects  which  could  lead  to  destabilization. 

To  recapitulate,  we  have  here  a  second-order  linear  oscillator  with 
stochastic  parametric  and  external  excitation.  We  can  write 

x  +  (an  +  a')x  +  (0  +  0')x  =  g 


where  aQ  =  <a>,  0Q  =  <0>,  a',  0',  and  g  are  statistically 

independent  stochastic  processes  with  given  statistics  to  second-order. 

Conventional  analyses  assume  white-noise  processes  which  is  neither 

2  2 

physical  nor  necessary.  Using  L  for  d/dt  ,  we  have 
Lx  +  OqX  +  0qX  +  a'x  +  0'x  =  g 
Let  R  =  aQ  d/dt  +  0Q  and  R  =  a'  d/dt  +0'  to  write 


M 


Wvlv 

m 


a 


m 


102. 


* 

» 


Lx  +  Rx  +  i?x  =  g 


or 


Lx  =  g  -  Rx  -  i?x 


which  is  immediately  solvable  by  decomposition  as  we  have  seen  from  previous 
examples.  All  that  is  necessary  is  operation  on  both  sides  with  L~\ 
decomposition  of  x  into  Zxn,  and  identification  of  the  Xq  term  to 
include  L_1g  and  the  initial  condition  terms.  It  remains  solvable 
even  if  we  add  a  nonlinear  term  Nx  (without  assuming  weak  nonlinearity). 


ACTIVE  DAMPING  OF  RESPONSE  OF  LARGE  SPACE  STRUCTURES 


Large  space  structures  are  necessarily  flexible  so  that  measures  must 
be  taken  to  minimize  structural  dynamic  response  to  orbital  corrections. 
The  problem  of  active  damping  to  transient  response  is  considered,  e.g., 
by  Chen  ( Int .  J.  Spacecraft,  21,  no.  5,  Sept-Oct  1984,  pp.  463-467).  Let 
us  consider  a  relevant  model  equation  derived  by  Chen: 

d^u/dt^  +  e|u|(du/dt)  +  u  =  0 
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are  not  commenting  on  this  particular  stiffness  control  scheme  but  only 

pointing  out  that  we  now  have  a  better  solution  procedure  for  any  such 

model.  This,  of  course,  was  the  objective  of  Phase  I.  In  further  work 

now,  we  can  implement  the  new  methodology  on  computers  for  application 

* 

to  real  structures  being  designed  by  aerospace  contractors. 


Also  we  have  seen  how  we  can  generalize  control  theory  to  apply  to 
nonlinear,  stochastic,  and  distributed  cases  by  analytic  rather  than 
numerical  procedures.  This  was  not  a  part  of  our  objective  but  it  turned 
out  an  unexpected  benefit  of  considerable  significance. 
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A  final  example  and  general  remarks: 

As  a  basis  of  discussion  and  pertinent  remarks,  let  us  look  at  a 
rather  general  nonlinear  partial  differential  equation  of  the  type  that 
could  arise  in  some  space  application.  Let  us  consider  the  (parabolic) 


equation 


ut  *  auxx  +  f<u’t.x)ux  =  0 


on  a  defined  finite  region  on  R  with  t  >  0.  Assume  a  is  a  constant, 
f  is  a  smooth  function  of  t,x,u,  and  the  initial/boundary  conditions 
are  given. 

Of  course,  to  obtain  a  quantitative  solution,  the  specific  form  of 
f  is  required  which  is  known  when  we  have  a  particular  model  to  solve.  Since 
each  model  will  be  different,  we  have  tried  to  show  how  any  model  can  be 
calculated.  When  f  is  given,  any  separable  terms  in  x  and  t  will  be 
designated  by  -g(x,t)  and  any  remaining  term  dependent  on  u,  and  multiplying 
u  can  be  written  N(u,uJ.  Let  l.  =  3/9t  and  LVy  *  32/3x2  and  write 

X  /N  V  AA 

Ltu  "  aLxxu  =  9  "  N<u’ux) 

The  decomposition  method  solves  for  each  linear  operator  term  in  turn; 


In  zero-gravity  environments,  heat  transfer  is  thought  to  occur 
only  by  radiation  or  conduction.  However,  it  is  possible  that  heat 
transfer  by  convection  can  also  occur,  e.g.,  in  heat  pipe  technology 
for  cooling  space  stations  sometimes  used  in  satellites.  This  can  orcur 
if  fluids  are  involved,  perhaps  in  machinery.  Rapid  heating  near  a 
boundary  can  cause  pressure  waves  producing  a  velocity  and  heat  transfer 
by  convection.  For  such  a  convection-diffusion  process,  we  get  a 
parabolic  equation  of  this  form. 
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Ltu  =  g  +  aLxxu  -  N(u,ux) 


Lxxu  =  "“'^9  +  a"lLtu  +  a"^N(u»ux) 

The  inverse  for  the  operator  l_t  is  an  integration;  similarly 

L  ^  is  a  two-fold  integration  since  L  represents  a  two-fold  differen- 

A  X 

tiation.  Since  lA.u  =  u  -  A  =  u  -  u(x,0)  and  L“!l  u  =  u  -  B  -  Cx, 
we  obtain 

u  =  A  +  L^g  +  aL^u  -  ^Vu.u,,) 

u  =  B  +  Cx  -  cr'l^g  ♦  <x-1L^Ltu  +  a'1!^  N(u,ux) 

A  is  the  initial  condition,  B  and  C  are  evaluated  from  the  remaining 
two  conditions.  We  add  the  two  equations  for  u  and  divide  by  two, 
obtaining  a  single  equation  for  u.  If  we  define 

uQ  =  (l/2){u(x,0)  +  B  +  Cx  +  L^g  - 

K  =  -(WHaCt\  +  a-’L^Lt) 

G  =  -(1/2HL"1  -  a_1L^} 

we  have 


u  =  Ug  +  Ku  +  G  •  N(u,ux) 

00 

The  solution  by  decomposition  is  u  =  £  un  with  the  given  Ug  and  the 
remaining  components  given  by 


u  , i  =  Ku  +  GA„ 
n+1  n  n 

for  n  ^  0,  where  the  An  are  the  polynomials  defined  by  the  author 
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in  [1]  for  N(u,utf),  i.e.,  N(u,uv)  =  T  An.  For  analytic  nonlinearities 
x  x  n=0  n 


for  which  the  A  can  be  determined  and  for  specified  initial/boundary 
''  00 


conditions  leading  to  a  convergent  series,  the  solution  is  u  =  £  u  , 


although  an  n-term  approximation  <|>n  =  J  u^  for  some  n  serves  as  a 


practical  solution,  usually  with  relatively  small  n,  as  seen  in  [1]. 


The  solution  can  also  be  made  by  operating  on  the  equation 


Ltu  "  aLxxu  =  9  '  N<u’ux)  with  (41  " 


Though  the  result  is  the  same  and  <j>t,  4>x  are  evaluated  from  the 
given  conditions,  writing  uQ  as 


uQ  =  (l/2)($t  +  4>x)  +  (1/2)(L'1  -  a-1L^)9 


means  the  result  is  not  limited  to  the  given  parabolic  equation  -  the 


derivatives  could  be  of  any  order  -  the  operator  L  represents  the  highest- 


order  derivative.  In  this  particular  case  we  used  the  double  subscript 


L  to  denote  a  second-order  derivative.  The  A  can  be  generated  for 
xa  n 


wide  classes  of  composite  nonlinear  functions  [1,2].  Thus  the  method 


is  quite  global 


When  solutions  are  obtained,  we  find  they  converge,  generally 


with  extreme  rapidity,  so  a  few  terms  of  u  =  £  up  suffice,  i.e., 

n-1 

we  use  a  practical  (n-term)  approximation  u  =  T  u . .  If  we  take  n 

£  1 


very  large,  this  approximation  approaches  £  u.  which  is  the  definition 

i=0  1 

of  u.  We  are  dealing  with  physical  systems  with  bounded  inputs,  analytic 


nonlinearities  for  which  the  An  can  be  determined,  and  specified 


initial/boundary  conditions,  assuming  integrability  of  the  forcing 
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function  -  all  physically  reasonable  conditions  to  get  convergence. 

In  stochastic  cases,  the  series  is  stochastic  and  from  the  approximate  u, 
we  can  determine  expectation  or  covariance. 
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INTRODUCTION  TO  THE  PROBLEM  OF  CONTROL 


The  control  problem  involves  not  just  solving  a  differential  equation 


-  the  state  space  equation  -  but  finding  a  solution  satisfying  a 


performance  functional  or  criterion  to  be  minimized  (such  as  least  fuel, 


least  co«t,  least  time,  etc.).  Bellman  has  discussed  the  differential 


equation  dy/dt  =  g(y,r),  y(0)  =  c  and  the  criterion  equation 


(T 

J(u,v)  =  h(y(t, ) ,  r(t,))dt,  +  k(y(T)) 
Jn  1  1  ! 


viewing  r(t)  as  a  point  in  function  space  and  setting  the  variation  of 


J  equal  to  zero  to  obtain  a  necessary  condition  for  the  determination  of 


r(t).  This  procedure  yields  the  Euler  equation  and  some  auxiliary 


conditions  which  convert  the  problem  of  minimizing  the  integral  functional 


to  that  of  solving  a  set  of  differential  equations  for  y  and  r.  These 


equations  are  nonlinear  and  subject  to  two-point  boundary  conditions. 


Since  the  decomposition  method  has  been  demonstrated  to  yield  solutions 


for  formerly  analytically  intractible  nonlinear  or  stochastic  equations. 


it  is  clear  that  a  viable  approach  to  the  general  control  problem  is 


possible.  The  difficulty  previously  was  not  the  formulation  of  the 


Hamiltonian  approach  but  the  inability  to  solve  stochastic  or  nonlinear 


equations  as  they  are  or  a  global  approach  to  partial  differential 


equations.  The  problem  has  now  been  shown  to  fit  in  the  format  of  the 


decomposition  method. 


Thus  one  can  consider  a  state  space  equation  of  the  form 


x  =  f(xp...,xn;  Uj,...,u  ;  t)  and  general  performance  criterion  functions 


J  extending  to  the  stochastic  and  multidimensional  cases  for  application 


to  space  structures.  The  next  section  considers  our  approach  to  the  analysis, 
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AN  OUTLINE  OF  OUR  APPROACH  TO  CONTROL 


Suppose  we  consider  a  nonlinear,  possibly  stochastic  which  we  want 
to  control  in  some  optimal  way.  For  a  linear  control  system  with  a 
quadratic  performance  index,  of  course  an  analytical  solution  can  be  made. 
Consider  the  state  equations 

x(t)  —  f  (x.j , . . .  ,X^ i  u i , . .  • , u^ ,  t ) 

i.e.,  a  set  of  n  nonlinear  differential  equations  with  x(t)  representing 
a  state  vector  with  n  components  f-|,...,f  ,  and  x(tg)  a  given  initial 
vector.  Define,  for  example  [4]  a  performance  functional  J(x,u,t)  given 
by 

rti 

J  =  4>[x ( t-j ) ,  t-|]  +  F(x,u,t)  dt 

t0 

where  <p  and  F  are  scalar  functions  with  necessary  smoothness  pr.perties. 
Let  p  =  [p-|,...,Pn]^  be  a  vector  of  Lagrange  multipliers  and  form  an 
augmented  functional 

J'  =  ^IXt,),  t,]  +  [F(x,u,t)  +  pT  (f-x)]  dt 

Ho 

Integration  by  parts  leads  to 

J '  =  4>  -  [pT  x]  j^1  +  [H  +■  pT  x]  dt 

to  H0 

with  H  defined  as 

H(x,u,t)  =  F(x,u,t)  +  pTf 

If  u  is  defined  on  tg  <  t  <  t^ ,  we  vary  u  and  find  the  variation 
6J 1  corresponding  to  6u,  leading  to  the  n  adjoint  equations. 


m 


p.  =  -  9H/3xi 

so  we  have  a  system  of  2n  nonlinear  differential  equations  with  two-point 
boundary  conditions.  Although  this  approach  has  been  discussed  by 
R.  E.  Bellman  and  many  others  perhaps  most  recently  in  [4],  analytical 
solution  has  usually  not  been  possible  except  by  numerical  methods.  We 
now  have  a  promising  and  potentially  valuable  alternative  since  such 
systems  of  nonlinear  differential  equations  have  been  solved  (even  for 
the  stochastic  and/or  multidimensional  cases)  in  an  analytic  approximation 
by  the  decomposition  method  [1-3]. 

Another  possibility  is  through  solution  by  decomposition  of  the  matrix 
Riccati  equation  which  appears  in  invariant  embedding  and  neutron  transport 
theory  as  well  as  modern  control  theory.  Consider 

R'(x)  =  B(x)  +  D(x)R(x)  +  R(x)D(x)  +  R(x)B(x)R(x) 

R(0)  =  0 

where  B,  D,  R  are  continuous  n  x  n  non-negative  matrices.  Suppressing 
the  argument  x,  we  have 

R'  =  B  +  DR  +  RD  +  RBR 

If  L  =  d/dx 

LR  =  B  +  HR  +  NR 

where  LR  =  R',  HR  =  DR  +  RD,  and  NR  represents  a  nonlinear  operator 
on  R.  Since  R(0)  =  0,  operation  with  L"^  on  both  sides  yields 

R  =  L”1 B  +  L-1HR  +  L_1NR. 
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Let  R  and  NR  be  written  in  terms  of  the  An  polynomials.  For  R  this 

oo  00 

is  equivalent  to  writing  R  =  £  R  .  For  NR  we  write  l  A  .  Identify 

n=0  n  n=0 

Rq  =  L_1 B ;  then 

R0  =  L-’B 

R1  =  L_1HR0  +  L_1A0 
R2  =  L~ 1  HR-j  +  L_1A1 


R  =  L-1HR  ,  +  L-1An  , 
n  n-1  n- 1 


for  n  >  1.  The  An  for  NR  are  given  by  [1] 


Aq  -  r0br0 


A-j  =  RqBR1  +  R^Rq 


A2  =  R-jBR-j  +  RqBR2  +  R2BRq 


A3  =  RQBR3  +  R3BRg  +  RlBR2  +  R2BRi 


a4  =  r2br2  +  r0br4  +  r4br0 

+  R^  BR^  +  R3BR1 


so  that 


V  L  B 

R1  =  L_1HR0  +  L-1R0BR0 


R~  =  L^HR,  +  L"1  (Ra.BR,  +  R,BRn) 


R3  =  L_1HR2  +  L"1(R1BR1  +  RqBR2  +  R^Rq) 


Finally,  since  HR  =  DR  +  RD 
R0  =  L-'B 


R,  -  L''(DR0  +  RqD)  +  L'1(R0BR0) 

R2  =  L‘1(DR1  +  R-|D)  +  L“1(R0BR1  +  R^Rq) 

R3  +  L-1(DR2  +  R2D)  +  L“1(R1BR1  +  R0BR2  +  R2BRq) 


i  *  « 

An  n-term  approximant  is  4>  =  £  R  which  approaches  R  =  £  R 

n  i=0  n  i=0  n 

as  n  ->  °°.  Thus  given  B,  D,  a  specific  R  can  be  calculated  to  a  desired 
approximation.  Accuracy  has  been  demonstrated  in  [5]. 
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FINAL  REMARKS 


When  we  study  a  physical  system,  we  must  simplify  the  actual  system 
with  an  idealized  model  represented  by  a  set  of  equations  which  can  then 
be  solved  by  available  techniques.  Thus,  modeling  is  always  a  compromise 
between  a  sufficiently  realistic  representation  and  tractability  of  the 
resulting  equations.  The  real  test  is  whether  the  resulting  solution 
yields  results  in  close  conformity  with  the  actual  physical  system.  As 
our  ability  to  solve  more  complicated  equations  increases,  we  become  less 
limited  in  making  our  models  physically  accurate.  Thus,  as  computers  have 
grown  in  capability  and  speed,  there  has  also  developed  an  increasing 
reliance  on  them.  This  dependence  has  in  some  cases  given  rise  to  an 
attitude  that  analytical  solutions  are  neither  possible  nor  really 
desirable.  It  is  indeed  correct  that  closed  form  solutions  are  generally 
not  possible,  and  that  everything  else  can  be  considered  "approximation." 
Yet,  all  modeling  is  approximation,  and  closed  form  solutions  mean  very 
little.  Linearization,  for  example,  is  a  general  procedure  and  has 
contributed  much,  but  the  linearized  model  is  not  the  same  problem  and 
may  lose  important  aspects  of  the  real  behavior.  When  we  consider  that 
closed  form  solutions  result  from  changing  the  physical  problem  until 
the  equations  are  sufficiently  simple,  we  realize  that  we  are  solving  a 
different  problem. 

We  have  demonstrated  that  the  decomposition  method  allows  physically 
realistic  and  accurate  solutions  of  the  types  of  vibration,  heating,  and 
control  problems  which  will  arise  in  the  analysis  of  large  space  structures 
planned  for  the  next  two  or  three  decades.  The  severe  requirements  imposed 
by  the  sizes,  weights,  necessary  stiffness,  and  specifications  on  thermal 
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and  mechanical  distortions  make  this  methodology  essential  because  of  its 
particular  features. 

The  chief  advantage  is  that  solution  of  nonlinear  problems  and/or 
random  problems  is  possible  without  linearization,  assumption  of  weak 
nonlinearity,  perturbation,  assumption  of  special  (and  unphysical 
stochastic)  processes,  or  classic  numerical  methods.  The  elimination  of 
the  restrictive  assumptions  and  limitations  of  usual  methods  will  reduce 
computational  effort  and  add  certainty  to  the  design  optimization  process 
of  rapidly  maneuvering  space  structures! 

It  is  essential  in  dealing  with  the  complex  interdependent  problems 
arising  in  contemplated  space  structures  that  we  make  physically  realistic 
solutions,  not  mathematical  solutions  of  models  which  have  been  made 
solvable  by  restrictive  assumptions.  Numerical  solutions  by  super¬ 
computers  have  brought  us  far  but  have  drawbacks  of  sometimes  excessive 
computer  time  and  a  lack  of  insight  into  functional  dependences.  We 
need  insight  into  performance  and  correct  solutions  of  the  physical  problems 
to  prevent  failure  in  operation.  The  decomposition  method  allows  us 
a  new  and  valuable  alternative  to  be  used  by  aerospace  contractors. 


So  far,  we  have  shown  in  the  analysis  of  problems  arising  in  space 
structures  that  the  decomposition  method  provides  efficiency  and  accuracy 
as  well  as  other  important  advantages.  These  advantages  are  insight  into 
functional  dependences,  physically  more  correct  and  realistic  solutions, 
and  a  considerable  potential  savings  of  computational  effort.  This  work 
applies  to  the  response  of  nonlinear  structures  involving  random  vibrations, 
fluctuations  of  parameters,  thermal  response  in  problems  involving  random 
or  temperature-dependent  thermal  conductivity,  radiative  transfer,  and 
control. 

Two  additional  avenues  have  become  apparent.  First,  a  major 
generalization  and  extension  of  control  theory  to  the  analytical  solution 
of  the  nonlinear  stochastic  multidimensional  case  has  now  become  possible. 

Second,  programming  of  the  decomposition  method  is  essential  now.  We 
believe  this  will  be  a  major  step  with  significant  advantages  over  numerical 
methods  now  in  use. 

Both  of  these  are  of  crucial  significance  to  programs  such  as  those 
contemplated  by  SDIO.  The  computer  implementation  of  decomposition  will 
allow  us  to  consider  problems  of  greater  complexity  and  general  contractor 
models.  The  control  theory  will  allow  us  to  do  the  problems  of  distributed 
nonlinear  stochastic  control  essential  for  space  structures.  Both  should 
be  pursued  as  separate  projects  to  be  brought  together  for  Phase  III  type  work. 

The  overall  program  is  intended  to  allow  the  solution  of  any  specific 


problems  of  interest  in  real  space  structure  models  of  interest  to  Air  Force, 
NASA,  or  space  station  contractors. 
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